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CALCULATION OP UNSTEADY TRANSONIC 
AERODYNAMICS FOR OSCILLATING WINGS 
WITH THICKNESS 
(COMPUTER PROGRAM) 


By S. Y. Ruo 

Lockheed-Georgia Company 


SUMMARY 

A computer program has been developed to account approximately for the 
effects of finite wing thickness in the transonic potential flow over an 
oscillating wing of finite span. The program is based on the original sonic- 
box program of Rodemich and Andrew, and accounts for the nonuniform flow 
caused by finite thickness by application of the local linearization concept. 
A brief description of each subroutine is given, and the method of input is 
shown in detail. A sample problem as well as a complete listing of the 
computer program are presented. 



INTRODUCTION 


This report is prepared as a supplement to NASA Contractor Report CR-2259 
entitled "Calculation of Unsteady Transonic Aerodynamics for Oscillating Wings 
with Thickness." The purpose of this supplement is to aid those who are 
interested in using the computer program. 

A major part of the computer program described herein is adapted from the 
original sonic-box computer program developed by Rodemich and Andrew..* In 
addition to its original capability, the program has been extended to include 
the thickness effect of a finite wing in unsteady sonic flow. This is 
accomplished in an approximate manner by employing the local linearization 
concept in conjunction with the known steady Mach number distribution or 
pressure distribution of the wing at its mean position. As the program is 
presently formulated, Mach number is assumed to depend only on the wing 
thickness; therefore, the local Mach number is the same at corresponding points 
on upper and lower wing surfaces. The flow is assumed to remain attached and 
shock-free. The program is limited to wings with unswept trailing edge and 
without control surfaces. The plane perpendicular to the wing planform and 
passing through the centerline chord is the plane of symmetry for both 
geometry and motion. 

The wing with thickness is converted into a wing without thickness in a 
transformed space in which the flow properties over the wing are computed. 

These transformed flow properties are then converted into the corresponding 
properties in the physical space and the generalized aerodynamic force 
coefficients are computed. 


* Rodemich, E. R. ; and Andrew, L. V.: Unsteady Aerodynamics for Advanced 

Configurations, Part II - A Transonic Box Method for Planar Lifting 
Surfaces. FDL-TDR-64-152, Part II, May 1965, Air Force Flight Dynamics 
Lab., Wright-Patterson Air Force Base, Ohio. 



SYMBOLS 


a 

nm 

AXY(I,J) 

AY(J) 


b 

d , , 
n'm' 

^(x.y) 

f(x,y) 

G-,(x,y), G 2 (y) 
i 
k 
L 

L. . 

10 

S 


coefficients in doublet potential polynomial, see equation (2) 
integral defined in equation (5-) 
integral defined in equation (5) 

area of ij-th box on the wing (i-th spanwise column and j-th 
chordwise row) 

reference length - wing centerline chord (dimension = L) 
coefficients in wing deflection polynomial, see equation (3) 
leading edge correction term, see equation (4) 
mode shape polynomial, see equation (3) 
leading edge correction terms, see equation (4) 

reduced frequency, cub/U^ 
unit of lengtn 

generalized aerodynamic force coefficient 
dimensionless wing planform area of full wing (reference 

,2s 

area = b ) 


T 

U 

00 

w (x,y) 

x,y 

x t (y) 

y 

‘'max 

(x,y) 

\Kx,y) 


unit of time 

freestream velocity (dimension = i/T) 
dimensionless downwash (reference velocity = U ) 
dimensionless cartesian coordinates (reference length = b) 
x-coordinate of wing leading edge 

dimensionless maximum semi-span (reference length = b) 

dimensionless velocity potential of doublet 

downwash at point (x,y) due to doublet of unit strength at 
origin, see equation (8) 

angular velocity (dimension = radian/T) 
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VARIABLES IN COMMON STATEMENT 


The variable names listed in the general COMMON statement in the MAIN 
program are described briefly in this section. The subroutine name in the 
parentheses following the description of a variable indicates where its value 
is generated or defined. The number appearing in the parentheses indicates 
the value of the variable for the corresponding condition described. In the 
following, N=1 and N=2 indicate, respectively, the real and the imaginary 
parts of a variable; i and n are related to the chbrdwise coordinate while 
j and m are related to the spanwise coordinate; * represents multiplication. 

A(N,j,i) influence coefficient, the integral in equation (6); the 

upwash at the center of a box caused by a unit doublet distri- 
bution over another box separated from the former by j number 
of boxes in spanwise direction and i number of boxes in 
chordwise direction, (P0T2) 

V(N,j,i) downwash at ij-th box, equation (?), (WVAL) 

S(N,j,i) velocity potential at ij-th box, (BOXPO) 

DA(k) input data, (DATRD) 

PS(N,n,m*M) inner summation of equation (l) including the integrals for 
M-th mode, (BOXP) 

DF(NEW,n,m*M) coefficients in wing deflection polynomial at M-th mode, (DRED) 
ML(NEW,i) number of boxes in i-th spanwise column, (SHAPE) 

AXY(n,m) integral defined in equation (5)> (SECT) 

AY(m) integral defined in equation ( 5 ), (SECT) 

XEDGl(k) x-coordinate of input leading edges of physical wing, (SHAPE) 

YEDGl(k) y-coordinate of input leading edges of physical wing, (SHAPE) 

IEDG(NEW) flag to identify the last section of the leading edge is 

parallel (l) or not parallel (0) to the freestream, (SHAPE) 

NS(NEW) number of leading edge segments, (SHAPE) 

EM(j,i) Mach number at the center of the ij-th box calculated from 

the fitted surface, (SHAPE) 

XEDG(k) x-coordinate of leading edges used in computation, (SHAPE) 

YEDG(k) y-coordinate of leading edges used in computation, (SHAPE) 
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/ 


HM(n,m) 

AREA 

D 

DI 

DH 

CK 

L 

M 

NEW 

IR 

IW 


coefficients in Mach number polynomial, (MRED) 
area of the physical, full wing planform, (SHAPE) 
length of box side, (MAIN) 

number (real) of boxes along centerline chord, (MAIN) 
one-half of the length of box side, (MAIN) 
reduced frequency 

number (integer) of boxes along centerline chord, (MAIN) 
deflection mode number, (MAIN) 

index for physical (1) or transformed (2) wings, (MAIN) 
computer read-unit number, (MAIN) 
computer write-unit number, (MAIN) 
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FUNCTION OF SUBROUTINES 


A brief description of the function of each subroutine is given in this 
section. Some equations are noted here for the purpose of identifying the 
steps in the sequence of obtaining the final results, i.e., the generalized 
aerodynamic force coefficients. It is suggested to refer to the report of 
Rodemich and Andrew* for the derivations of these equations. The freestream 
flow is in the positive x direction and the apex of wing is at the origin. 
Spanwise direction is designated by y. Since the plane passing through the 
centerline chord and perpendicular to the wing planform is the plane of 
symmetry of wing geometry and motion, the input of wing geometry, deflection, 
and steady-flow Mach number distribution needs to be made only for the portion 
of the wing where x > 0 and y > 0. It is assumed that the camber, twist, and 
mean angle of attack of the wing are all small; and the difference in local 
Mach number at the corresponding points of the upper and lower wing surfaces 
can be neglected. As the computer program is presently formulated, the 
motion is limited to the symmetric mode. This limitation is imposed by the 
expression used in the present program to evaluate the influence coefficient 
of upwash due to a unit doublet distribution. 


In the following, the variable name used in the computer program to 
represent a certain quantity mentioned in the description of a subroutine is 
indicated in the parentheses following such quantity. 


MAIN - A controlling routine; also performs the final calculation of the 
generalized aerodynamic force coefficient which is written for full wing as 


L ij S , d n'm' 

° n' ,m' 


£ / 

n,m L 

+ik^ )y 2 ( m+m )p(x,y)dxdy 

S 

-n'-| //> n -1 ) y 2 ( m+m ^F(x,y)dxdyJ , 


2(m+m,) F(l,y)dy 


( 1 ) 


wnere S = area of the wing planform (full wing) 
k = reduced frequency, 

dn'm' a nm are res P ec ^ ve ^-y coefficients of the wing 

deflection and velocity potential polynomials, i.e. 


*see footnote on page 2 
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( 2 ) 


c P ) (x,y) = I § a nm xU y 2m ? (x»y), n,m = 0 — 4 

f(x,y) = £, S, d n , m , x n y 201 , n',m' =0^4 (3) 

The function i!‘(x,y) in equation (2) is defined as 

*’(x,y) = & 1 (x,y) . G 2 (y) (4) 


where the expressions for G^ and G^, depending solely on the shape of the 
leading- edge, are 


G-,(x,y) = 


My) = 


V 


2 2/ v 

x -x^ (y), 


if y = 0 at x = 0 


or 


y x-x^(y) , if y / 0 at x = 0 


ffc) 


or 


, if ^ = 0 at i « 1 
dx 


. if § / 0 at x " 1 


These relationships are depicted in figure 1, and x (y) is the leading edge 
of tne wing (see fig. 1). 

DATKD - data input routine; for a new wing, the data array (BA(k)) is 
cleared every time hut for the same wing, the old value is used if a new 
value is not entered (blank datum will not affect the old value). 


SHAPE - wing geometry routine; approximates the wing planform with a 
grid of square boxes and also calculates the number of boxes (ML(REW,i)) 
along each spanwise column (i-th) on both physical and transformed wings; 
computes the area of the full wing planform (AREA) and the surface-fitted 
(hM(n,m)) Mach number (EM(j,i)) at the center of each box (ij-th) on the 
pnysical wing; increases the number of leading edge segments (XEDG(i), 

YKDG(i)) by subdividing the input leading edge segments (XEDGl(i), YEDGl(i)) 
for the physical wing so that the transformed wing can be better approximated; 
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performs the coordinate transformation of the wing with the aid of the known 
Mach number distribution (EM(j,i)) on the physical wing and also checks the 
possible creation of an artificial wake in the transformed wing resulting 
from an uneven distortion of the wing tip. 

MACH - Mach number subroutine; computes the mean local Mach number at a 
given point (x,y) on the physical wing from the known polynomial (HM(n,m)) 
which was fitted to the input data for Mach number distribution; the Mach 
number at corresponding points on upper and lower surfaces is assumed to be 
the same. 

FORCI - controlling routine for the calculation of integrals in the 
expression for generalized aerodynamic force coefficient (Eq. (l)); obtains 
the values of 


AXY(I,J) = 1 JJ x^ 1-1 )y 2 ^ J "’ 1 ^P(x,y)dxdy 
S 

AY(J) - \ f 


( 5 ) 


y2( J_1 ) F ( 1f y) dy> 


X=1 


where the integrals are evaluated in subroutine SECT. Since these integrals 
depend only on the wing geometry, their values are calculated once and used 
for all frequencies and modes both with and/or without thickness. 

SECT - carries out the calculation of the integrals in equation ( 5 ) 
using the six-point Gaussian quadrature 

n 


n l n n 

g(x,y) dxdy = EE h.h.g(x. ,y .) 
. _ i-i >1 ■> J 


and 


jf 


n 


g(x)dx = £ h.g(x. ), 
i=1 1 1 


in which n = 6 is assumed. 

P0T2 - evaluates the integral in the expression relating the downwash 
and the velocity potential 


2 cp 0 / y.- T l)d|dTi = w(x ,y ) 

L 9 J 1 o g 

i'j* 


( 6 ) 
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where <¥> 


o . , . , 


= magnitude of velocity potential at the center of box IL, 


w(x.,y.) = downwash at the center of box B. . 
v i ij 

”(lf + ikf ) at (x i’ y P 


( 7 ) 


^(x . -£,y .-p) = downwash at (x.,y.) due to doublet of unit 
strength at 


ik 1 ( 1.. r .s ( y r , > ) il 

■ 2 ^ -—2 2 ik L ( vhn^rjrJi ' 




( 6 ) 


the integral in equation (6) is integrated around the edge of each box , 
centered at (x^i,yji) once for every frequency and the result (A(lJ,j,i)) 

is used in the calculation of the velocity potential (S(M,j,i)) in both 
physical and transformed spaces for all modes under consideration at the 
same frequency. 


Clh - calculates sine and cosine integrals. 


DEED - fits the wing deflection for each mode (M) of motion in both 
pnysical and transformed spaces with a polynomial of the form shown in 
equation (3) expressed in terms of the dimensionless coordinates and stores 
the coefficients (DF(i'IEW,n' ,m'*M)); if the data of wing deflection are given 
at a number of points on either the physical or the transformed wing, these 
daca are fitted into the desired polynomial by the least-square surface 
fitting method; for the physical wing, the wing deflection in each mode can 
be given as input either in the form of the desired polynomial coefficients or 
deflection at a number of selected points; for the transformed wing, the wing 
deflection in a particular mode is obtained from a number of selected points 
on the physical wing of known deflection (DF(l ,n' ,m'*M)), the deflection at 
a point on the transformed wing is unaltered from that of the corresponding 
point on the physical wing but its spanwise coordinate in the transformed 
space is modified by a factor of the Mach number from that of the physical 
wing. 

WVAL - calculates the complex (real-N=1, imaginary-N=2) downwash 
( w (N,j,i)) at the center of each box (ij-th) for a particular mode (M) and 
reduced frequency (CK) using the coefficients of wing deflection polynomial 
(DF(i\EW,n' ,m'*M)) for the physical or the transformed wing. 

BOXPO - solves a set of simultaneous equations, equation (6), relating 
the complex (real-N=1, imaginary-H =2 ) velocity potential (S(N,j,i)) and the 
known complex downwash (W(b , J , i ) ) of a particular mode (ii) at the center of 
each box (ij-th) to obtain the velocity potential of the physical or the 
transformed wing. 
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BOXP - fits the complex (real-N=1, imaginary-N=2) velocity potential 
distribution (S(N,j,i)) calculated in BOXPO for physical or transformed wing 
with a pre-determined polynomial surface (C(k,N)), as shown in equation (2), 
for both real and imaginary parts; for the case with thickness effect, the 
complex velocity potential on the physical wing is obtained from the fitted 
complex velocity potential (C(k,N)) for the transformed wing by selecting a 
number of points on the physical wing, computing the complex velocity 
potential of the corresponding points on the transformed wing with the known 
polynomial and transforming the potential to that of those points on the 
physical wing, then fitting this indirectly obtained potential with the pre- 
determined polynomial surface (C(k,N)); also calculates the values of the 
inner summation (PS(N,n,m*M)) of equation (l); prints out the downwash 
(W(H , j , i ) ) calculated from the deflection polynomial (DF(NEW,n,m*M)), the 
computed potential (S(N,j,i)), the fitted potential (PSI(N,j)) or the lifting 
pressure (PR(W,j)) calculated from the fitted potential, etc., at the box 
centers, if the controls (DA(1003) to DA(1010)) are flagged. 

MRED - fits the input steady-state Mach number distribution or pressure 
coefficient distribution on the physical wing with a polynomial surface of 
the form shown in equation ( 3 ) expressed in terms of the dimensionless 
coordinates for Mach number; the input for Mach number can be either in the 
form of the desired polynomial coefficients or the Mach number or the pressure 
coefficient at a number of selected points on the physical wing; for the 
latter case, a least-square surface-fitting method is used to obtain the 
desired polynomial (HM(n',m')). 

LSQUA - method of the least-square surface-fitting routine. 

MSIMER - finds solutions of the simultaneous real number equations. 

MSIMEC - finds solutions of the simultaneous complex number equations. 
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FLOW CHART 


The flow chart of the computer program is presented in figure 3* The 
function of those controlling variables appeared in the flow chart is as 
follows. 


variable 

VALUE 

FUNCTION 


1 

Indicates case without thickness effect 

NEW 


2 

Indicates case with thickness effect 


0 

Indicates first frequency for a wing 

DA(26) 


1 

Indicates additional frequency for the same 
wing 


0 

Calculates cases with and without thickness 
effects 

DA(44) 

1 

Calculates case without thickness effects 
only 


2 

Calculates case with thickness effects only 

M 

< DA(28) 

Counter of the modes of deflection 

DA(28) 

< 10 

Total number of modes to be considered 

-lWD 


Index of the storage location in data array 
DA 
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DA(26) = 



15 


Figure - Flow chart 
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At tne start of execution, as it calls DATED, the data array, DA is 
cleared by a "+" sign entered in the first column of tne first data card. 

The wing geometry, reference velocity, first frequency, data for the first 
mode of wing motion and the optional set of controls (Da( 1001 )-DA(l014)) are 
entered with a sign at the first column of the last card of this data 
set. 


The physical wing planform is approximated by a grid of square boxes in 
SilaPh and those integrals depending only on the wing geometry are evaluated 
in f'ORCl and SECT. Those integrals involving the frequency and wing geometry 
are calculated in P0T2. The wing deflection of this particular mode is 
rearranged into a pre-determined polynomial surface in DEED. Subroutine 
LSviUA is not called when the wing deflection is entered as the desired poly- 
nomial coefficients. With the known downwash, computed in WVAL with the aid 
of the polynomial for deflection, the velocity potential at the center of each 
box is calculated in BQXPG. This potential is fitted into a polynomial in 
BOX.?. The unsteady lifting pressure may be calculated and printed out from 
BOX? if desired. This ends the calculation of one mode of wing motion at one 
frequency. If there are more modes to be treated, a new set of wing deflec- 
tions for each mode is read in separately each time and the corresponding 
velocity potential is computed until all modes are completed. Then the 
generalized aerodynamic force coefficients are calculated and printed out. 

This completes the computation for the case without thickness effect. 

If the case with tnickness effect of the same wing at the same frequency 
is desired, the steady local Mach number data are read in, either in a poly- 
nomial form or the actual values (steady pressure coefficient may be used) at 
a set of selected points on the wing, and they are fitted with a pre- 
determined polynomial surface in MEED. In subroutine SHAPE the Mach number 
at the center of each box over the pnysical wing and at the end points of 
eacn leading edge segment is computed using MACH with the aid of the poly- 
nomial for Mach number and the wing planform is then transformed; also the 
number of the straight line segments representing the transformed wing leading 
edge is increased so that the distorted wing can be better approximated. 

Since the frequency and the original wing have not changed, subroutine EORCI 
and P0T2 are skipped. There is no need to read in the wing deflection again. 
Each mode is treated as before except the wing is now considered to be in the 
transformed space. Tne wing deflection is redistributed on the transformed 
wing. This is accomplished by computing the wing deflection at a set of 
selected points on the physical wing, transforming the coordinates of those 
selected points and then fitting the wing deflection with a polynomial in 
terras of the transformed coordinates in DEED. The velocity potential com- 
puted in BOXPO for the transformed wing is converted and surface-fitted into 
that of the physical wing in BOXP as each mode is treated. The calculation 
of generalized aerodynamic force coefficients is performed in the physical 
space. This completes the calculation for a wing oscillating at different 
modes at one frequency without and with thickness effect. 

The whole process is repeated for each additional frequency except the 
subroutines SHAPE, F0KC1, DEED and 'MEED are skipped. A (DA(44)) flag nay be 
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used to control the calculation for the case either with or without thickness 
effect or for both. 

There is no limit to the numbers of different frequencies or different 
wings that can be handled in a single execution. A blank card placed behind 
the last card of the data sets will automatically terminate the execution. 
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iLiPUT GUIDE 


Data input is made through the subroutine DATED using the one dimension- 
al array DA with a size of 101 5. The allowable maximum number for some of the 
input data as indicated below may be changed if the dimension of the corres- 
ponding storage array and computational operations are also changed accord- 
ingly. The layout of the array DA(k) as it is presently used is as follows: 

1-12: Title 

13-22: Mode title 

23 : Frequency, (cycle/sec) 

24: Centerline chord length, (ft or meter) 

25: Speed of sound of the freestream, (ft/sec or meter/sec) 

26: (0) - indicates the frequency is the first one for a new wing 

(1) - indicates the frequency is the additional one for the 
same wing 

2?: Humber of boxes along centerline chord (maximum 50 ) 

23: Number of deflection modes (maximum 10) 

2i'. Humber of sections of leading edge to be given (maximum 6) 

30-42: Coordinates of points on the leading edge, (ft or meter) 

(in the sequence of y Q , , y. } , x 2> y 2 » ... , x^, y^) 

43 : Humber of boxes allowed for upstream influence (if this 

location is left blank or assigned a zero, it will assume 
DA(43)=DA(27) and in no case DA(45) DA(’27) is allowed). 

44: (0) - indicates to calculate cases with and without thickness 

effect 

12 

45: Indicator to suppress calculation of potential for a mode 

(0) - no suppression 

(1) - suppression 


- indicates to calculate case without thickness effect only 

- indicates to calculate case with thickness effect only 


46 - 7 C: Coefficients of the deflection polynomial (in the sequence of 

d U0’ d 10’ d 40’ d 01’ d ir ’** ’ d 41» •'* ’ d 04’ d 14’ *'* ’ d 44 ; 
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yi-95: Coefficients of the Mach number distribution polynomial (in the 

sequence of ^.qq, ^-jq» •••* ^q» ^oi* d n’ **** ^ 41 * *^ 04 * 

d 14» d 44) 

96: Indicator of the type of wing thickness effect input 

(1) - pressure coefficient 

(2) - Mach number 

97: Number of points at which pressure coefficient or Mach 

number to be given 

98 : Number of points on which deflections to be given 

99: 1 plus the maximum power of x in surface fit of deflection 

(maximum 5 ) 

2 

100: 1 plus the maximum power of y in surface fit of deflection 

(maximum 5) 

101 -700 : Deflection data for a maximum of 150 points (in the sequence 

of x, y, deflection and weighting factor) 

701-1000: Pressure coefficient or Mach number data for a maximum of 100 

points (in the sequence of x, y and pressure coefficient or 
Mach number) 

The remaining part of DA array is used for the control of intermediate 
results print out. When the latter is desired, a non-zero number should be 
entered at a location in the DA array corresponding to the information wanted. 

1001: Coefficients of surface-fitted deflection polynomial (no 

thickness effect) 

1002: Coefficients of surface-fitted deflection polynomial (with 

thickness effect) 

1003: Value of upwash calculated (no thickness effect) 

1004: Value of upwash calculated (with thickness effect) 

1005: Value of potential calculated (no thickness effect) 

1006: Value of potential calculated (with thickness effect) 

1007: Coefficients of surface-fitted potential polynomial (no 

thickness effect) 

1008: Coefficients of surface-fitted potential polynomial (with 

thickness effect) 
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1009: Surface-fitted values of potential and lifting pressure (no 

thickness effect) 

1010: Surface-fitted values of potential and lifting pressure (with 

thickness effect) 

U11: Coefficients of surface-fitted Mach number distribution 

polynomial 

1012: Surface-fitted value of Mach number distribution 

1013: Coordinates of sections of the leading edge of the transformed 

wing 

1014: Value of deflection on transformed wing 

1015: Not used 


The format of the input data card is (A1 , A5, 16, 11A6, A2). The first 
field is for the control of clearing the data array, DA, for a new wing (+) 
and the control to indicate the end of the set of data (-). The second field 
is tne indicator for the type of daxa, either numeric (blank) or alphameric 
(AnPHA). The third field is the designator for the relative location in the 
data array of the first numuer to follow in the fourth field. If this field 
is left blank, or a zero is entered, the execution will be terminated. The 
fourth field is for five consecutive input data each occupying 12 columns 
plus 0 blame columns at the end. All the fixed point numbers are right- 
adjusted and the decimal point for the floating point number must be included. 
If an input datum is left blank, no change at the storage location for that 
particular datum in the data array will occur unless the set of the input 
data is for a new wing. 
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SAMPLE CASE 


A typical data cards set-up and output for an aspect ratio 1.5 delta 
wing having an elliptic lateral cross-section with lOfe thickness ratio 
performing pitch about its apex and plunge are presented. 


Input 

The input format is (A1 , A5, 16, 11A6, A2). 


1 2 3. 4 5 6 7 

12 34 5678 90 123456 7890 1234567890 12 34 56 7890 1234 5 6 7890 12 34 56 7890 1234 56 7890 12 IDENTIFICATION 


♦ALPHA 

1ASHECT 

RATIO 1.5 

delta wing 

(TAU=0.10> 

1 


23 

1.591549407 

10.0 

1000.0 

3 

ALPHA 

13PLUNGE 




2 


27 


40 

2 

1 

4 


30 


0.0 

10.0 

3.75 

5 


H6 


1.0 



6 

ALPHA 

13PITCH 

AriOUT HOOT 

leading edge x=o.d 

7 


46 


0.0 

0.1 


8 


96 


1 

26 


9 


7U1 


0.1 


0.140232 

10 


704 


0.3 


0.135125 

11 


707 


0.5 


0.132175 

12 


7iU 


0.7 


0.129943 

13 


713 


0.9 


0.128084 

14 


716 


1.3 


0.124980 

15 


719 


1.7 


0.122340 

16 

„ 

722 


2.1 


0.119963 

17 


725 


2.5 


0.117746 

18 


72b 


2.9 


0.115627 

19 


731 


3.3 


0.113564 

20 


734 


3.9 


0.110510 

21 


737 


4.5 


0.107428 

22 


740 


5.1 


0.104243 

23 


743 


5.7 


0.100872 

24 


746 


6.3 


0.097215 

25 


749 


q.9 


0,093129 

26 


752 


7.5 


0.088395 

27 


755 


7.9 


0.084697 

28 


750 


6.3 


0.080342 

29 


7ol 


8.7 


0.-074970 

30 


764 


9.1 


0.067825 

31 


767 


9.3 


0.063062 

32 


770 


9.5 


0.056823 

33 


773 


9.7 


0.047661 

34 

- 

776 


9.9 


0.029651 

35 

- 

23 

3.103098014 


1 

36 

- 

23 

6.366197620 


1 

37 


38 


Card 1 : title of the case under consideration. 

Card 2: title of the first mode of deflection. 

Card 3: first frequency (cycle/sec), centerline chord length (ft), 

reference velocity (ft/sec).' 
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Card 4 : number of boxes along the centerline chord, number of deflection 

inodes, number of total leading edge segments of the wing. 

Card 5s spanwise coordinate (ft) of the first section of the leading edge, 
chordwise and spanwise coordinates (ft) of the next section (the 
sequence is y Q , x^ , y^ — e.g. , see figure 2 ). 

Card 6 : first mode of deflection f = 1.0, 

the " - " sign indicates the end of the group of data cards to be 
read at this stage. 

Card 7: title of the second mode. 

Card 8 : second mode of deflection f = O.lx. 

Card j : identification of the type of input regarding the wing thickness 

effect (pressure coefficient for this case), number of points on 
the wing this information to be given. 


Cards 10 to 35: 

chordwise and spanwise coordinates (ft) of a point on the wing, and 
the pressure coefficient* at this point. 

the " - " sign on the last card indicates the end of the group of 
data cards to be read at this stage. 


Cards 3° and 37s 

additional frequencies for the same wing, one card is read in at 
one time. 


Card 38 : blank card to make an exit from the computer. 


*Alksne, A. Y. j and Spreiter, J. R. : Theoretical Pressure Distributions 

on Wings of Finite Span at Zero Incidence for Mach Number Near 1. 
d ASA TR R- 88 , i 960 . 
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Output 


aspect ratio i.s delta wing (tau-o.ioi 

.. MO BOXES ALONG ROOT CHORD - .. . ROOT CHORO LENGTH « ln.CC FT - 

REOUCEO FREQUENCY- .100 FREE STREAM VELOCITY .■400**00 FT/SEC 

. - FREQUENCY ■ 1. 592 + 00 CYCLE/SEC - 

- MODE NO. I PLUNGE - — 

MODE N0»— 2 PITCH-AB0UT-. ROOT-LEAO INt-EOGE *»3.0 


GENERALIZED FORCES (NO THICKNESS EFFECT) 

MOOES 

_PRE5» OEFL« PEAL PART IMAG .PART. ABS.. VALUE — pHASeaNGLE 

I l -6.97765-PS -2. 90975-01 2.90975-01 -9C.0IS9 

1 2 - 1 .95306-09 -1.58646-01 1.58666-01 _ -90.0706 

2 1 -2.90*59 + 00 _ : --209303*81 2.92199+00 -179.3289.. - 

2 2 -1.59065+00 -1.77008-01 1.60097+00 .173.6502 


MODE. NO* 1 - - : - - -i - — . 

MOOE NO. 2 - - 

GENERALIZED FORCES (WITH THICKNESS effect) 

MODES 

PRES. OEFL. . .._ . REAL PART IMAG PART APS. -VALUE . pHAS£ ANGLE 

1 1 1 .559 11-03 -2.6029 2-01 2 .5029 7-01 -89.693| 

1 2 1.22693-03 -1.63389-01 1 *63399-01 -89.56*9 

2 1 -2.99939 + 00 -2.46717-01 2.51363 + 00 .173.9087 

2 2 s4 .6JC35+UO- -4-. 987.99=04 1.6 9 2 1 2+00 -47-3.0968 - 
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aspect ratio us delta wing (TAU»0.lO) 

SO BOXES ALONG ROOT chord .. ROOT CHORD LENGTH . . • In.CP FT 

REOUCEO frequency ...» .2CC - — . ..... FREE -STREAK VELOCITY-* lOQji.OC FT/SEC- 

. . FREQUENCY • 3 • 1 83 + 00 CYCLE/SEC - - 

MODE N 0 • - I . . . .... . . 

. - MODE NO. 2 . . . .. 

GENERALIZED FORCES (NO THICKNESS EFFECT ) 

MOOES 

PRES* -OEFL* *, REAL .PART JMAG..PART ABS». VALUE pHASe ANGLE.. 

I I S.33I9C-03 -N. 77153-01 S. 77173-01 -B9.N79a 

1 2 3.912C1-03 -3.|S|66-3i 3. I 9189-01 -B9.3778 

2 ... | - -2.399S8 + C0 -H.95SS3-CI 2.«5.O23*00 - -1*8.3390 

2 2 - - 1 .58 1 75 + 00 -3. *9772-01 I. *2939*00 -I**,8s2o 


HOPE NO. 1 -- - 

MODE NO. 2 - - 

GENERALIZED FORCES (WITH THICKNESS EFFECT) 

MODES 

PRES. OEFL. .. REAL PART IMAG PART APS* VALUE PHASc ANGLE 

- 1 1 — *.30539-03- -S.9B9ia-CI 9. 9 89 58-0 3 -89.2759 - 

12- s. 77907-03 -3.25*37-01 - 3.25*72-01 -89,1592 

2 I -2.S05S5-00 -1 - ~ -5.29752-01 2, 55595*00 -.1*8.0380 

2 2 -1 .*319 1*00 "3.937 33-C 1 U*7«73*nO. - _.l*i.H35S- .... 
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-- ASPECT ratio 1»S OELTA « I NG ( T AU«0 • 1C 1 1 

.... 90 BOXES ALONG ROOT CHORD-- ROOT CHORD LENGTH 

REDUCED FREQUENCY- » .900 FREE. -STREAM-VELOCITY 


FREQUENCY ■ 6. 346 + 00 CYClE/SEC 

MODE NO. 1 

-MODE NO.- 2- - 

GENERALIZED FORCES (NO THICKNESS EFFECT) 


MOOES 

-PRES.-OEFL. REAL PART IM*G -PART ABS». YALUE.. . 

1 1 3.C8793-02 -9.98383-01 9.«8BB5-0I 

12 2.96913-02 -6.29898-01 6.25336-01 

- 2.~~ -1 2 * 909 95 + 00 — ... •1.00687+00 —2. 60675+00 — 

2 2 -1.58739+00 -7.53332-01 1 . 75703*00 


MODE NO* l 
MODE NO. 2 


GENERALIZED FORCES (WITH THICKNESS EFFECT) 

MODES 

PRES». DEFL. REAL PART .IMAG PART AES . VALUE 

— l 4 3 + 33909-02 - -9.87578-01 9.88132*01 

1 2 2.97565-02 -4.93576-01 6.99052-01 

— 2 - 1 -2.99756 + 00 -l. 05458+00 2.71186+00 

Z 2 1.63073+00 — n-7 .83179-01- 1 .80909+03 


Ip.OC FT 

- tono.oc FT/SEC - 


PHASE ANGLE 

- - 88,1351 

-8Z.737J 

. 157.2784 

- 159 . 6)15 


PHASE ANGLE 

—-88,0909 

-87.7971 
- .157,06*9 
.159.3968 
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COMPUTER PROGRAM LISTING 


C MAIN PROGRAM 

COMMON A (2*100,50) ,W(2»50*50) »0A(1015> *PS(2» 5,50 ) »DF (2, 5, 50 ) 

1 , ML(2»50) » Axy (9, 9) » AY (9) * XEDGI (8 ) * YEOgX ( g) 

2 ,IE0G(2) ,NS(2) *EM(50*50) ,XE06<32) *YE0«C32)*HM(5*5) 

3 , AREA, D*DI*DH*CK*L,M, NEW, IR#XW 
DIMENSION G(10) 

DATA Z/1H / 

IR=5 

IW=6 

WRITE(IW,55) 

C READ DATA FOR THE ACTUAL WING 

100 CALL DATRO(DA) 

NEW=1 

TESTl=OA ( 1001 )+OA( 10o3) +DA( 1005) +DA( 1007) +DA(1009) 

CK=0A(23)*0A<24)/0A(?5)*6.283l8531 

DI=DA(27) 

L=OI 

0=1. 0/DI 
DH=0.5*D 
WRITE 

0 (IW, 60) (DA(I) ,1=1,12) 

110 IF (L) 600,600,120 
120 IF (50-L) 600,130,130 
130 WRITE 

0 ( IW ,65) L,DA(24)*CK,DA(25>*DA<23> 

IF (DA (26) ) 160,150,160 
150 CALL SHAPE 

IF(NEW.EQ.2) GO TO 180 
CALL FORCI 
160 LIM=ML(NEW,L) 

IF (LIM-50) 170,170,650 
170 LIM2=2*LIM 

NFLNS = DA (43) 

IF < NFLNS. E«.0> NFlNS = L 
LPOT = MIN0(L, NFLNS) 

CALL POT2(100*LIM2,LPOT,CK,0) 

180 CONTINUE 
M=0 

K=DA(28) 

GO TO 230 

C PRELIMINARY CALCULATIONS ARE FINISHED, 

C THE NEXT SECTION IS GONE THROUGH FOR EACH MODE. 

200 IF (DA (26) ) 230,210*230 
210 IF (NEW.EG.2) GO TO 230 
CALL OATRD(DA) 
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230 K=K-1 
M=M+1 

IF(TESTl.LT.l.O) SO TO 250 

o"* ( i w , 1 o > m REPRODUCIBILITY OF THE 

10 format ( mi lsx , shmode no. 13) ORIGINAL PAGE IS POOR 

GO TO 270 
250 WRITE 

0 ( IW» 15 ) M 

15 FORMAT (1H0.15X.8HMODE NO.*l3) 

270 CONTINUE 

WRITE 

0 (IW.16) (DA(I> .1= 13.22) 

16 FORMAT (1H+.30X.12A6) 

IF (OA(26)) 310.290*310 
290 CALL DRED 
300 G(M)=DA(45) 

310 IF (G(M) ) 380.320.380 

320 IF(IFIX(DA(44)).EG.2.AND.NEw.EQ.1) GO TO 360 
CALL WVAL 

IF CML(NEW.l)) 330.370.330 
330 LIM2=2*ML(NEW*1) 

DO 350 J=1.LIM2 

350 W(J.1.1)=W(J. 1.1 >*2.0/3.14159265 
C LEADING EDGE CORRECTION 
370 CONTINUE 
CALL BOXP 

380 IF (K) 400.400.390 
390 IF (M-10) 200.400.40Q 

C FINAL SECTION OF PROGRAM - COMPUTATION OF GENERALIZED FORCES 
400 CONTINUE 

IF(IFIX<DA(44) ) .EQ.2. AND.NEW.EQ.l) GO TO 480 
IF (TEST1.LT. 1.0) GO TO 410 
WRITE( IW.55) 

WRITE 

0 ( IW.60 ) (DA ( I ) ,1=1,12) 

WRITE 

0 (IW.65) L,DA(24) ,CK,0A(25) ,DA(23) 

410 WRITE (IW.20) 

IF (NEW .EQ. 1 ) WRITE ( Iw , 22 ) 

IFINEW.EQ.2) WRITE ( Iw » 2*+) 

WRITE( IW.25) 

20 FORMAT (/1H010X.18HGENERALIZED FORCES) 

22 FORMAT (1H+.28X.22H (NO THICKNESS EFFECT)) 

24 FORMAT (1H+.28X.24H (WITH THICKNESS EFFECT)) 

25 FORMAT (1H0.5X.5HM0DES/4X, 11HPRES. DEFL. ,8X,9HREAL PART, 10x»9HIMAG 
SPART.10X.10HABS. VALUE , 6X7 1 lHPHASE ANGLE) 

AC=8. 0/AREA 
DO 470 Mlsl.M 
IF (G(M1) ) 470.430.470 
430 DO 450 M2=1.M 
S1=0.0 
S2=0.0 
N1=5*(M1-1) 
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N2=5*(M2-1) 

DO 440 J=l*5 
J1=J+N1 
J2=J*N2 
DO 440 1=1*5 

Sl=Sl*PS ( 1 » I * Jl ) *DF ( 1 , 1 * J2 ) 

440 S2=S2+PS(2»I*J1)*0F(1#I#J2) 

S1=AC*S1 

S2=AC*S2 

S3- SQRT(S1**2+S2**2) 

S4=57 . 29578*ATAN2 { S2 . SI ) 

450 WRITE 

0 <IW*30)M1*M2.S1,S2*S3*S4 

30 FORMAT (1H02I6* 1P3E19. 5* 0P1F16.4 ) 

470 CONTINUE 
480 WRITE ( IW*55) 

IF(DA(26) .6T.0.) 60 TO 520 
DO 500 1=13*22 

500 DA ( I )=Z 

C CHECK IF TRANSFORMED wING HAS BEEN CALCULATED 

C IF NOT* 60 TO SHaPE TO TRANSFORM THE WXN« 

520 IF(NEW.EQ.2.0R.IFlX(0A<44) ) .EQ.l) SO TO 100 
NEW=2 

IF (DA(26) ) 180*530*160 
530 CALL MRED 
GO TO 150 
C ERROR EXITS 
600 IPR=27 
610 WRITE 

0 <IW*35)IPR 

35 FORMAT(1H010X»8HBAO 0ATAI4) 

GO TO 700 

650 WRITE <IW*40) 

40 FORMAT (1H010X*42HLATERAL LIMIT ON NUMBER OF BOxK EXCEEDED* ) 
700 STOP 
55 FORMAT (IH1) 

60 FORMAT (1H010X*12A6) 

65 FORMAT ( 1H0*10X* I2*23H BOXES ALONS ROOT CHORD* 

1 15X.22HROOT CHORD LENGTH =*F8.2*3H FT/ 

2 1H0*10X»19HREDUCED FREQUENCY =*F6.3* 

3 15X*22hFREE STREAM VELOCITY =»F8. 2* 7« FT/SEC/ 

4 1H0*10X»11HFREQUENCY =* lPEll «3* 10H CyCLE/SCC) 

END 
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SUBROUTINE DATRO(OATA) 

c card-read subroutine «datrd(Data(i))* 

DIMENSION DRBU{12) *DATA(1) #DDRBU(lO) 

DATA ATEST/5HALPHA/*DTEST/lH /,ETEST/1H-/,STEST/1H*/ 
DATA Z/1H / 

IR=5 
IW=6 

1 READ ( IR * 2) EMIN. ALP * IND, <DRBU< I ) , Isl , 12) 

2 FORMAT (Al*A5*l6*llA6*A2) 

IF ( IND.EQ, 0 ) SO TO 20 
IF (EMIN.NE.STEST) SO TO 105 
NEW WINS IF COLUMN 1 CONTAINS A PLUS SISN 
INITIALIZATION OF DATA ARRAY 
DO 101 1=1 * 1015 

101 DATA ( I ) =0 • 0 
DO 102 1=1*22 

102 DATA ( I ) =Z 
DO 103 1=104,700,4 

103 DATA ( I ) =1 .0 
105 CONTINUE 

IF (ALP.EO.ATEST) 60 TO 9 
IF (ALP.NE.DTEST) GO TO B 
C NUMERIC CARD 

DO 3 1=1*10 
DDRBU ( I ) =DRBU ( I ) 

3 CONTINUE 

DECODE (DDRBU* 990) (ORbU ( I ) * 1=1 * 5) 

DO 5 1=1*5 
IF(DRBUd) )4*6*4 
C TEST FOR BLANK FIELD 

6 IF(SIGN (1.0*DRBU(I) )) 5*5*4 

4 DATA ( IND)=DRBU( I ) 

5 IND=INDU 
60T0 11 

C ALPHA CARD 
9 DO 10 1=1*10 

DATA ( I NO ) =0RBU ( I ) 

10 INO=INO +1 

11 IF (EMIN.NE.ETEST) SO TO 1 

C RETURN IF COLUMN 1 CONTAINS A MINUS SIGN 
13 RETURN 

C END OF DATA CAROS 
20 WRITE( IW*995) 

STOP 

C BAD CARD 
8 CONTINUE 
WRITE 

0 ( IW*993) EMIN* ALP* IND. IDRBU( I ) * 1=1 * 12) 

WRITE (IW*991) 

STOP 

990 FORMAT (5E12*0) 

991 FORMAT (38H BAD DATA ON THIS CARD. JOB TERMINATED > 
993 FORMAT (12H0CARD IMAGE=A1 * A5, 16, 12A6) 

995 FORMAT (///lH *10X,17HNO MORE DATA CARD/) 

END 



subroutine shape 

COMMON A (2* 100 #50) #W (2 *50*50 ) *DA (1015) »PS (2 #5# 50 ) »DF (2# 5# 50 ) 

1 *ML{2#50)#AXY<9#9)#AY(9)»XED6I(8)*YEDGl(8> 

2 #IEDG(2) #NS (2) *EM(50#50> *XEDG(32) *YEDg<32> #HM<5»5> 

3 *AREA,D»0I#DH*CK*L#M#NEW»IR*IW 

c following common is a temporary storage in this Routine 

COMMON ER(5> 

ERRR=l,E-06 

IEDG(NEW)=0 

IF(NEW,EQ.2) GO TO 200 
NS<1)=DA(29) 

IF (NS(1) ) 750 # 750 #100 
100 IF(NS(l)-6) 105# 105*750 
105 NSP=NS(1)+1 

IF (DA(24) ) 755*755# 110 
110 DO 115 I=1*NSP 

XEDG(I)=DA(2*I+27)/DA(24) 

115 YEDG(I)=DA(2*I+28)/DA(24) 

XEDG(1)=0.0 

XEDG (NSP*1> =1.0 

YEDG(NSP+1>=YEDG(NSP) 

DO 120 1=1*8 
XE0GI(I)=XEDG(I) 

120 YEOGI(I)=YEDG(I) 

GO TO 500 

C FOR MODIFIED WING 

C COMPUTE MACH DISTRIBUTION ON ORIGINAL WING 

200 DUM=0. 

DO 220 1=1 »L 
X=DH+FLOAT(I-l)*D 
JML=ML ( 1 # I ) 

DO 220 J=1*JML 

Y=DH+FLOAT(J-l)*D 

CALL MACH(X#Y*EMLOC*DUM*HM*l) 

220 EM(J#I)=EMLOC 
K=1 

NSP=l,+DAC27)/2. + l.E-07 
D2=2»*D 

C DEFINE EDGES OF MODIFIED WING 
DO 240 I=2*NSP 
XEDG ( I ) =FLOAT ( 1-1 ) *D2 
IF(XEDG( I ) .GT.XEDGI (K) ) K=K+1 

YEDG ( I ) =YEDG I ( K-l > ♦ ( YEDG I (K> -YEDG I C K-l ) ) * ( XEDG ( I ) -XEDGI ( K-l ) ) / 
$ (XEDGI (K) -XEDGI (K-l) ) 

240 CONTINUE 

IF(IFIX(DA(27))-2*(NSP-1) I 710*255*250 
250 NSP=NSP+1 
NS12=NS( 1 ) +2 
XE06 ( NSP ) =XEDGI ( NS12 ) 

YEDG ( NSP > = YEDG I ( NS 1 2 ) 

255 CONTINUE 
Kl=2 
K2=NSP 

NSl2=NS(lJ+2 
DO 300 I=2*NS12 
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OO 260 K=K1iK2 
K3— K 

IF(ABS(X£DGI(I)-XEDG(K> ) .LT.ERRR) GO TO 290 
IF(XEDG(K) . GT.xEDGI ( I ) ) GO TO 270 
260 CONTINUE 
60 TO 290 
270 K4=K2+K3 

DO 280 J=K3'K2 
K=K4-J 

XE06(K+1)=XEDG(K) 

280 YEDG(K*1)=Y£0G(K) 

NSP=NSPn 

XEDG(K3>=XEOGl(I> 

YEDG(K3)=YEDGI(I) 

K2-K2+1 
290 K1=K1+1 
300 CONTINUE 

IF (DA ( 1013) ) 310*320,310 
310 Kl-1 

NS(NEW)=NSP 

NS12=NSP+1 

WRITE 

0 (IW*40) Kl* NEW, NS (NEW) , NEW » IEDG (NEW ) 

WRITE 

0 (IW»50) (K*XEDG(K) , yEDG(K) ,K=1»NS12) 

320 CONTINUE 

IF(NSP.GT.32) GO TO 720 

C CHECK IF MODIFIED WING HAS A FOLD-OVER EDGE 
C TRANSFORM Y-COORDINATE 
DO 350 IsifNSP 
YMP=0. 

Y1=0, 

Y=OH 

330 CALL MACH(XEDG(I) * Y *EMLOC«DMDY » HM* 2) 

YM=Y*EMLOC 

ymd=ym-ymp 

IF(Y1.6T.Y> ymd=-ymd 
IFUMD.GE.O.) GO TO 335 
ER(1)=XEDG(I) 

ER<2)=Y 

ER(3)sEML0C 

ER(4)sDMDY 

ER(5)=YMP 

IERsS 

C ACCEPT A FOLD-OVER OF NOT MORE THAN 1/5 OF A BOX 
IF(ABS(YMD).GT.0.4*DH) GO TO 7l5 
IERROR=330 
WRITE 

0 ( IWf 20) IERROR* (ER( J) » J=l* IER) 

335 IF(Y.EQ.YEDGd) ) GO TO 340 
Y1=Y 
Y=Y+0 
YMP=YM 

IF(Y.GT.YEDGd) ) Y=YEDG(I> 

GO TO 330 

340 YE06d)=YEDG(I)*£ML0C 
350 CONTINUE 

IF(DA(1013) ) 360 » 370 *360 
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360 Kl=2 
WRITE 

0 (IWt40) Kit NEW# NS (NEW) #NEW» I EDS (NEW) 

WRITE 

0 { IW»50) (KtXEDG(K) rYEDG(K) »K=1*NS12) 

370 READJUST THE POINT DISTRIBUTION TO DEFINE NS(2) 

K1=0 

1CH£CK=0 

TMP=0. 

DO 450 I=2#NSP 
YMP=AMAX1 ( YMP » YEDG ( I ) ) 

IF ( YEDG (I> -YEDG ( 1-1 > ) 400 # 420#440 rner 

ACCEPT SMALL WAKE BEHIND TRANSFORMED WINS LEADlNS ED6E 
OF WING# IF NOT MORE THAN 1/5 OF A BOX 
400 ER(1)=XEDG(I) 

ER(2)=YEDG(I> 

ER(3)=XEDG(I-1) 

ER(4)=YEDG(I-1) 

ER(5)=YMP 

1ER=5 _ , _ 

IF<YMP-YEDG(I).GT.0.4*DH) GO TO 705 
IERROR=400 

write 

o ( IW»20) IERROR,(ER(J)#J=l*IER) 

YEDG(I) = YMP 

420 IFU.EQ.ICHECK+1) GO TO 430 
lCHECKsl 
GO TO 440 
430 lCHECKsl 
K1SK1+1 
440 K=I-Kl 

XEDG(K)=XEOG(I) 

YEDG(K)=YEDG(I) 

450 CONTINUE 

1F(ABS(1.-XEDG(K)) .GT.ERRR> 60 TO 700 
IF ( YEDG ( K ) -YEDG ( K-l ) ) 705*460t470 
460 NS(2)=K-2 
GO TO 480 
470 NS(2)=K-1 
480 NSP=NS(2)*1 

IF (DA ( 1013) ) 490*495,490 
490 Kl=3 
WRITE 

0 (IW»40) KltNEWtNS(NEW) *N£W* IEDG(NEW) 

WRITE 

0 ( iWt 50) (KtXEDG(K)tYEOG(K) »K=1»NS12) 

495 CONTINUE 
YISYEDG(I) 

500 CONTINUE 

IF (Yl) 760*520»520 
520 K=0 
X=DH 

IF(NEW.EQ.2) GO TO 530 
AREAsO.O 

530 DO 600 I=1*L 
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540 IF <X-XEDG(K+1>> 550,550*570 
550 ML(NE*»I>=0.54DI*(YED&<K)4F*<X-XEDG(K) )/G) 
GO TO 600 
570 K=K+1 

G=XEDG<K+1)-XEDG(K> 

F=YEDG(K+1>-YED6(K> 

IF (G) 765,580,580 
580 IF <F) 770,590,590 
590 IF(NEW.EQ.2> SO TO 540 

AREA=AREA+G*(YEDG(K+1)+YEDG{K) ) 

GO TO 540 
600 X=X+D 

DO 610 1=2, L 
K=I 

IF(ML<NEW»I>.LT.MLlNEW,I-l>> GO to 725 
610 CONTINUE 
K=L 

IF (ML (NEW, K) «GT*50) GO TO 725 
IF(NEW«EQ.1> GO TO 650 
IF(OAU012)> 620,650,620 
C PRINT OUT MACH DISTRIBUTION 

620 WRITE( IW»60) 

DO 640 1=1, L 
JL=ML(1,I) 

IF(JL.EQ.O) GO TO 640 
WRITE 

0 ( IW,70> I 

JLP=JL/6 

IF(JL-6*JLP.NE.0) JLP=JLP+1 
DO 630 J=1 ,<JLP 
630 WRITE 

0 ( IW, 80) ((J1,EM(Ji,I)),U1=J,UL,ULP) 

640 CONTINUE 
650 CONTINUE 

NS12=NS(NEW) 

IF(XED6(NSP).GT.1.) GO TO 765 
IF ( XEDG ( NSP ) -1 . 0+ERRR > 680 , 670 , 670 
670 IF (YED6(NSP)-YEDG(NS12>> 680,680,690 
680 IEOG(NEW)=l 
690 CONTINUE 
RETURN 

700 IERROR=700 
ERCl)sXEDG(K) 

ER(2)=YEDG(K) 

IER=2 
GO TO 740 
705 IERROR=705 
GO TO 740 
710 IERR0R=710 
ER(1)=0A(27) 

ER ( 2 ) =NSP 
IER=2 
GO TO 740 
715 IERRORS715 
GO TO 740 
720 IERR0R=720 

ER ( 1 ) =XEDG ( NSP ) 



ER(2)=YEDG(NSP) 

£R<3)=NSP 
IER=3 
GO TO 740 
725 IERR0R=725 
ER(1)=N£W 
ER(2)=K 

ER(3)=ML(NEW.K) 

ER(4) =ML(NEW*K-1 ) 

IER=4 
740 WRITE 

0 (1**20) IERROR, (PR(I)*I=1*IER) 

STOP 

750 IPR=29 

GO TO 790 
755 I PR-24 

GO TO 790 
760 IPR=30 

GO TO 790 

765 K= MlN0(K*NS(l) ) 

IPR=2*K+29 
GO TO 790 
770 IPR=2*K+30 
790 WRITE 

0 ( IW* 10 ) IPR 

10 PORMAT(1H0*10X.17hSHAPE — BAD DATA* 15) 

20 FORMAT ( 1H0* 10X*38hBAD NUMBER IN SHAPE NEAR STATEMENT N0.»l5* 

1 /1H »15X,1P5E14.6) 

40 FORMAT (/1HO*5X*3HNO., 12 *43H REDISTRIBUTION OF WIN® LEADING EDGES* 

1 NS(,I1*4H) = *12. 8H, IEDS(*U»4H) = *11) 

50 FORMAT ( 4(6X* 13* 1P2E11 >4) ) 

60 FORMAT ( 1H1 * 5X.30HLOCAL MACH NUMBER DISTRIBUTION/) 

70 FORMAT (1H0.5X* 12. 6HTH ROW) 

80 FORMAT (1H *5X»6(2x* I3.1PEI3.5) ) 

STOP 

END 


SUBROUTINE MACH(X» Y.EMLOC *DMOY*HM*K ) 
C CALCULATE LOCAL MaCH NUMBER 
DIMENSION HM(5*5) 

Y2=Y*Y 

YP=1.0 

EMLOCSO.O 

DMDY=0.0 

DO 20 J=l*5 

XYPsYP 

CsFLOAT(J) 

DO 10 1=1*5 

EMLOC=EMLOC+HM ( I * J) *XYP 
IF(J.EQ.5.0R.K.EQ.l) GO TO 10 
DMDY=DMOY+C4HM ( I * J+l ) *XYP 
10 XYP=X*XYP 
20 YP=Y2*YP 

DMDY=2.*Y*OMDY 

RETURN 

END 
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SUBROUTINE FORCI 

C CONTROL SUBROUTINE FOR CALCULATION OF INTEGRALS UstO 

C TO FINO GENERALIZED FORCES „ _ ^ 4 

COMMON A(2»100*50) ,W(2*50»50> »0A(l0l5) »PS(2*5*50) »DF(2»5»50) 

1 ,ML(2»50) »AXY(9»9)/AYl9)»XEDGI(8)»YEDGl(ft) 

2 * I£DG(2) »NS(2> »EM(50»50> #XEDG(32> »YE0G<32) #HM{5*5> 

3 , area, o«oI»oh»c*ul.m. new. ir*iw 

C FOLLOWING COMMON COMMUNICATES BETWEEN «F0RCI* AND ’SECT’ 

COMMON YMAX2.N 
NSS=NS(1)+1 

YMAX2=YEDG ( NSS ) WYEDG { NSS ) 

NSS=NS(1) 

DO 3 J=l»9 
DO 2 1 = 1 » 9 

2 AXY 1 1 » J) =0 • 0 

3 AY(J>=0.0 

IF (OA(30J) 4»5»4 

4 N=0 

C SECT DOES THE CALCULATIONS FOR EACH SECTION OF THE PLANFOrM 

5 DO 7 I=1»NSS 

IF (Y£OG(l)-YEDG(I+l) ) 6»7»7 

6 Ns I 

CALL SECT 

7 CONTINUE 
RETURN 
END 


SUBROUTINE SECT 


axy(i»j> is the Integral over the planform of 

X**( l-l)*Y**(2J-2)*P<X» Y)*Q(Y) 


IN TERMS OF THE EQUATION OF THE LEAOING EDGE X = X0(Y) 
P(X»YJ = SQRTF (X**L-X0**L) 

WHERE L = 1 OR 2 

Q( Y) : 1 OR SQRTF ( l-( Y/YMAX )** 2 ) 

DEPENDING ON THE WING SHAPE 


AY(J) IS THE INTEGRAL OVER THE TRAILING EDGE OF 
Y**(2U-2)*P<1»Y)*Q(Y) 

COMMON A (2* 100 » 50) ,W(2*50»50> » DA (1015 1 »PS(2,5»50) » OF (2* 5 *50 1 

1 ,MU2f50)»AXY(9»9).AY(9)»XEDGI(8)fYEDGng) 

2 , IEDG(2) » NS ( 2 ) ,EM(50»50) ,XEDG(32) >YED6<32> »HM(5>5> 

3 ,AREA,D»DI »OH»CK»L, M»NEW, IR. IW 

FOLLOWING COMMON COMMUNICATES BETWEEN iFORCl« AND 'SECT* 

COMMON YMAX2»N 

DIMENSION U(6)#H(6) . . 

COMPUTES THE CONTRIBUTION TO AXY, AY OF THE SECTION OF WING 
BOUNDED BY A SEGMENT OF THE LEADING EDGE# THE TRAILING EDGE. 
AND TWO LINES ON WHICH Y IS CONSTANT 

A SIX POINT GAUSSIAN FORMULA IS USED FOR THE INTEGRATION 
OVER EACH VARIABLE 
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U(1)=0. 61930959 
U (2) =0.83060469 
U<3>=0. 96623476 
U(4 ) =0.03376524 
U(5)=0. 16939531 
U(6)=0, 38069041 
H(1)=0. 23395697 
H(2)=0. 18038079 
H< 3) =0.08566225 
H ( 4 ) =H ( 3 ) 

H(5)=H(2) 

H(6) =H ( 1 ) 

GAUSSIAN POINTS AND WEIGHTS FOR THE INTERVAL <0>1) 
X2=XE0G (N+l ) 

X1=XEDG(N) 

Y2=YEDG(N+1) 

Y1=YEDG(N) 

IF (N) 4.2.4 
2 X1=0.0 
Y1=0.0 
4 DY=Y2-Y1 
DO 19 J=1.6 
V=U(J) 

G=H( J) *DY 

IF (Y2**2-YMAX2) 7.6.7 
b G=2.0*V*G 
V=V*V 

7 Y=Y2-V*DY 
X0=X24V*(X1-X2) 

X0P=1.0-X0 
G=G* SQRT(XOP) 

YQ=Y*Y 

IF ( IEDG(NEW) ) 8.9.8 

8 G=G* SQRTU.0-YQ/YMAX2) 

9 E=2.0*X0P*G 

IF (DA (30) ) 10.10.11 

10 G=G* SORT (1.0+X0) 

11 DO 17 1=1.6 
U2=U(I)**2 
X=X0+X0P*U2 
F=E*H( I)*U2 

IF (DA (30 ) ) 12.12.13 

12 F=F* SQRT ( X+XO ) 

13 YP=1.0 

DO 16 M=1 . 9 
XP=YP 

DO 15 IL=1.9 

AXY(IL»M)=AXY(IL.M)+XP*F 

15 XP=X*XP 

16 YP=YQ*YP 

17 CONTINUE 
YP=1.0 

DO 18 M=1 .9 
AY(N)=AY(M)4YP*G 

18 YP=YQ*YP 

19 CONTINUE 
RETURN 
ENO 
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SUBROUTINE POT2(M2»Mo,NO#CK,D) 

THE VELOCITY FIELD OF A UNIFORM DOUBLET DISTRIBUTION 
OVER A BOX IS COMPUTED AT ALL POINTS AT WHICH IT wlLL BE 
NEEDED AND STORED IN THE ARRAY A IN COMMON 

MO# NO CONTROL THE NUMBER OF VALUES COMPUTED 

M2 IS THE RANGE OF THE SECOND SUBSCRIPT IN THE ARRAY# 
DIMENSIONED A(2»M2,N2), BUT TREATED HERE AS AN ARRAY 
WITH TWO SUBSCRIPTS 

DIMENSION A (2# 1 ) 

COMMON A 

M-MO 

N=N0 

DK=CK*D 

DK2=DK**2 

M1=M-1 

DK8S0K2/8.0 

DK4=2.0*DK8 

DK12=DK2/12«0 

CM=0.5 

DH=DK*0.5 

DM=0.5*DH 

D0=2.0*DK 

DDM=DD 

D1=0.25*DK2 

B5=DK2/24.0 

DO 3 1=1 #M 

01 = 0.0 

B4=2.0/DM 

B2=B5/B4-DH 

83=-0.5*B5 

D3=DH*B4+B5 

D4=DK8*B4 

D04S2.0*D4 

CN=I.O 

K=I 

C3=0.0 

C4=0.0 

C7=0.0 

C8S0.0 

DO 2 J=1 * N 

A1=DM/CN 

Cl=CM* COS(Al) 

C2=-CM* SIN( A1 ) 

C5=CM*CIN(A1,C6> 

C6=-CM*C6 

C9=Cl-C3 

ClO=C2-C4 

C11=C5-C7 

C12SC6-C8 

A(l»K)=B3*C9-B4*ClO-B5*C3-Bl*Cll-B2*Cl2 
A(2»K ) =B4*C9+B3*CiO-B5*C4+-B2*Cll-Bl*Cl2 
23 C3=Ci 
C4=C2 
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C7=C5 

C8=C6 

Bl=Bl-01 

B3=83-D3 

B4=B4-04 

04=04+004 

CN=CN+2.0 

2 K=K+M2 
CM=CM+1#0 
OMsOM+DDM 

3 DDM=DDM+DD 
DO 5 IL=1 • 2 
Kl = l 

DO 5 J=i»N 
DO 4 1 = 1 » Mi 
K=K1+M-I 

4 A(IL.K)=A(IL,K)-A(IL.K-1> 
A(lLfM)=2,0*A(IL,Kl) 

5 K1=K1+M2 
CM=0.0 
DM=0.0 
DOM=OK 

00 12 1=1 r M 

C7=0.0 

C8=0.0 

C9=0.Q 

C10S0.0 

P1=0.0 

P2=0«0 

CN=1 « 0 

B6=0.5*DK12 

K=I 

00 10 J=1»N 

A1=CH/CN 

A2=0M/CN 

IF (A1-0.2J 7*7»8 

7 Bl=2.0-Al**2/3,0 
B2=-DK/(6.0*CN> 

SO TO 9 

8 83= SIN(A1)/Ai 
B1=2.0*B3 

B2=<B3- C0S(A1))/A2-0H/CN*B3 

9 B3= C0S(A2J/CN 
B4= SIN(A2)/CN 
C3=B1*B3+B2*B4 
C4=B2*B3»B1*B4 
B5=0H*CN 
C1=B5*C4-2.0*C3 
C2=-2.0*C4-B5*C3 
C5=C1-C7 
C6=C2-C8 
P3=P2-B6*CN 

P4=P3+2 . 0*OK12* (CN-l . 0 > 

A(i#K)=A(l#K) +C5-P1*C6+P3*C3-P4*C9 

A<2»K)=A(2#K>+C6+Pl*C5+P3*C4-P4*ClO 

P1-P1+0H 

P2=P2+CN*DK4 

CN=CN+2.0 

C7=C1 
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C8=C2 
C9=C3 
C10=C4 
B6=B6+DK12 
10 K=K+M2 
CM=CM+DK 
DM=DM+DDM 

12 0DM=DDM*DD 

D3=CK/ (2. 0*3. 14159265) 

M1=M2-M 

K=1 

A1=0,0 
DO 14 J=1»N 
C1=D3* SIN(Al) 

C2=-D3* COS(Al) 

DO 13 1=1, M 

DFE =A(1,K)*C1+A(2,K)*C2 
A(2,K)=A(2,K)*C1-A(1,K)*C2 
A(1»K)=DFE 

13 K=K+1 
K=K+M1 

14 Al-Al+DH 
RETURN 
END 


SUBROUTINE DRED 

COMMON A(2»100.50) ,W(2»50»50> »DA(10l5) ,Ps( 2, 5, 50 > ,DF (2,5, 5n ) 

1 , ML(2,50),AXT(9,9),AY(9) , XEDGI (fl) , TED© I (fi) 

2 • IEDG12) » NS ( 2 ) ,EM(50»50) ,XEDG(32) , YEDg(32) ,HM<5»5> 

3 ,AREA,D,0I»DH,CK«L,M,NEW,IR»IW 

C FOLLOWING COMMON COMMUNICATES BETWEEN *DRED* AND »LSUQA* 

COMMON XQ ( 150 ) , YQ C 150 ) *DEF8 ( 150 > , WTQ ( 150 ) , CoE ( 5 • 5) 

NCONTSNEW+IOOO 

IF (NEW. EQ. 2) 60 TO 400 

NP=DA(98) 

IF (NP) 730,170,100 

C A POLYNOMIAL FOR THE DEFLECTION IS FITTED To VALUES 

C OF DEFLECTION AT GIVEN POINTS. 

100 NX=DA(99) 

NY=DA(100) 

IF (NX) 710,710,105 
105 IF (NY) 720,720,110 
110 IF (150-NP) 730,115,115 
115 IF (NP-NX) 710,120,120 
120 IF (NP-NY) 720,125,125 
125 CONTINUE 
KP=1Q0 

DO 140 IPsl,NP 
XQ(IP)=DA(KP*1)/DA(24) 

YQ( IP) cOA(KP+2) /DA (24) 

DEFQ(IP)=DA(KP*3) 

WTQ(IP)=DA(KP+4) 

IF (WTQ(IP) ) 740,740,140 
140 KP=KP*4 

C USE LEAST SQUARE METHOD TO CURVE-FIT DATA 
IF(DA(NCONT) ) 150,160,150 
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150 WRITE 

0 ( IW» 30) M 

WRITE ( I Wf 10) 

160 CONTINUE 

CALL LSQUA ( NP , NX , NY » NCONT ) 

GO TO 200 
170 YP=1.0 
K=1 

DO 190 J=l*5 
XYP=YP 

DO 180 I=l»5 
C0E(I»U)=XYP*DA(K+45> 

K=K+1 

1 8 0 XYP=XYP*DA(24> 

190 YP=YP*DA(24>**2 
200 K=25*<M-1) 

DO 220 1=1 ,25 
K=K+1 

220 DF(l,K,l)sCOE(I.l) 

RETURN 

- FOR MODIFIEO WING 

1 - CALCULATE DEFLECTION AT 100 OR LESS ON THE ORIGINAL WING 

2- TRANSFORM X* Y, DEF TO THE TRANSFORMED PLANE 

3- CURVE-FIT THE DATA 
400 MM=5*(M-1> 

J1=MM+1 

J2=MM+5 

NTOT=0 

DO 410 I=1»L 
410 NTOT=NTOT+ML(l,I) 

IF (DA ( 1014) ) 420*430,420 
420 WRITE ( IW, 50 ) 

430 CONTINUE 

NJ=FLOAT(NTOT)/100. ♦ 0.5 

IF(NJ.LT.l) NJ=1 

NPT=0 

NP=0 

K=0 

DO 520 I=l»NTOT*Nj 
NPSNP+1 

440 IF(I.LE.NPT) GO To 450 
K=K+1 
NPTP=NPT 
NPT=NPT+ML(1,K) 

GO TO 440 

450 XQ(NP)=DH-*FLOAT(K-l)*D 

jy=i-nptp 

YO(NP)=DH+FLOAT(JY-l)*D 

Y2=YG(NP)*YQ(NP) 

DEFQ(NP)=0. 

WTQ(NP)=1.0 

YP=1.0 

DO 470 J=J1'J2 
XYP=YP 

00 460 IX=1*5 

DEFQ { NP ) =DEFQ ( NP > +DF ( l , I X * J ) * X YP 
460 XYP=XYP*Xtt(NP) 

470 YP=YP*Y2 

IF (DA (1014) ) 480,490,460 
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480 Y2=YQ(NP> 

YP=DEFO(NP> 

490 CONTINUE 

YQ(NP)=YQ(NP)*EM(JY,K) 

C DEFQ(NP>=DEFQ(NP)/(Em(JY*K)+1.E-10) 

DEFQ(NP)=DEFQ(NP) 

IF (DA (1014) ) 500*520,500 
500 WRITE 

0 ( IW,60) I »K» JY »NP, XQ(NR) » YQ (NP) * Y2 »DEFQ ( NP) »YP,EM(JY*K) 

520 CONTINUE 

C USE LEAST SQUARE METHOD TO CURVE-FIT DATA 

IF(DA(NCONT) ) 550,560,550 
550 WRITE 

0 ( IW *40) M 

WRITE l IWi 10) 

560 CONTINUE 
C 

CALL LSQUA ( NP , NP * NP * NCONT ) 

00 580 J=l*5 
DO 580 1=1 '5 

580 DF ( 2 * I * MM+ J ) rCOE ( I * J ) 

RETURN 
710 IPR=99 

GO TO 750 
720 IPR=100 
GO TO 750 
730 IPR=98 

GO TO 750 
740 IPR=K+4 
750 WRITE 

0 ( IW * 20 ) IPR 

STOP 

10 FORMAT (lH010x*56HcOMpUTED DEFLECTION = SUM OF DEF(N,M)*X**(N-i)*Y» 
1* ( 2M-2 ) /1H010X* 54H ( IN DIMENSIONLESS COORDINATES - DISTANCE/CHORD L 
2ENGTH ) /1H09X . lHN7x * 1HM16X*8HDEF (N»M) ) 

20 FORMAT (1H0*10X,16hDRED — BAD DATA* 15) 

30 FORMAT (lH0*8X*26HpHYSICAL PLANE — MODE NO. ,13) 

40 FORMAT (IHO. 8X »29HTRAnsF0RMED PLANE — MODE NO. ,13) 

50 FORMAT (1H1, 5X.42HTRANSFORMEO OEFLECTlON INFORMATION — DrEO/) 

60 FORMAT ( 1H , 5X*4l5*lP8E12.5) 

ENO 


SUBROUTINE LSQUA ( NP * NX * NY tNR I TE ) 

C CURVE FITS WITH LEAST SQUARE METHOD 

COMMON A (2* 100, 50) ,W (2*50*50) *DA(1015> *PS(2*5*50) »DF (2,5*50) 

1 , ML (2* 50 )» Ax Y (9*9), AY (9) »XEDGI (8 ) *YEDgI(8> 

2 * IEDG(2) ,NS(2) *EM(50*50) ,XEDG(32) *YEDg<32) ,HM(5*5> 

3 ,area,d*di*dh*ck»l,m,new,ir*iw 

C FOLLOWING COMMON ONLY COMMUNICATES WITH THE CALLING ROUTINE 
COMMON XQ<150> ,YQ(15O)tOEFQ(150> *WTQ(150) ,C0E(5*5) 

C FOLLOWING COMMON IS USED AS A TEMPORARY STORAGE IN THIS ROUTINE 
COMMON IXB(5)*B(25) *C(25,25) *G(25) 

RITE=DA(NRITE) 

MX= MINO (NX , 5) 

MY= MINO (NY ,5) 
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IY=MY 
IXYrMX+MY 
00 2 J=1 * 5 
00 1 1=1 >5 

1 COE ( I . J ) =0 . 0 

2 IXB(J)=MX 
NC=MX*MY 

3 00 5 1=1 >NC 
DO 4 J=1»NC 

4 C ( I * J)=0.0 

5 BlI) =0.0 

DO 11 1P=1 » NP 
X=XQ (IP) 

Y2=YQ(IP)*YQ(IP> 

DEF=DEFQ( IP) 

WT=WTQ (IP) 

YP=1.C 

K=1 

DO 8 J=l»IY 

XYP=YP 

JX=IXB(J) 

DO 7 1=1. JX 

G(K)=XYP 

XYP=X*XYP 

7 K=K+1 

8 YP=Y2*YP 

DO 10 1=1 »NC 
DO 9 J=1»NC 

9 C(I»J)=C(DJ)+G(I)*G(J)*WT 

10 6(I)=B(I)+G(I)*DEF*WT 

11 continue 

K=MSIMER(25.NC.1»C»B) 

IF (K-l >22.22.15 

15 DO 16 1=1. I Y 

ip=im-i 

IF (IXB(IP)+IP-IXY) 16.17.17 

16 CONTINUE 

ixr=ixr-i 

GO TO 15 

17 IXB(IP)=IXB(IP)-1 

IF { 1XB( IP) ) 18. 18.19 

18 IY=IP-l 

19 NC=0 

DO 20 1=1. IY 

20 NC=NC+IXB(I) 

GO TO 3 

22 K=1 

DO 23 J=1 . IY 
JX=IXB( J) 

DO 23 1=1. JX 
C0E(I.J)=BIK) 

23 K=K*l 

IF (RITE) 61.66.61 
61 CONTINUE 

DO 64 J=l» IY 
JX=lXB( J) 

00 63 1=1. JX 
WRITE 

0 


( IW.42) I. J»C0E( I. J) 



REPRODUCIBILITY OF THE 
ORIGINAL PAGE IS POOR 


63 CONTINUE 

64 CONTINUE 
66 CONTINUE 

42 FORMAT (3X.2l8»lPE25»5) 
RETURN 
END 


C 


C 


C 

c 


SUBROUTINE WVAL 

EVALUATION OF THE UPWASH ARRAY „ 

COMMON A <2» 100.50 ) .W<2»50. 50) »DA< 1015) *PS<2»5»50) »DF(2»5»50> 

1 ,ML(2.50).AXY<9.9).AY(9>»XEDGl(8)*YEDGl(8> 

2 #IE0G(2).NS<2) .EM{50.50)»XEDG<32)»YEDg<32>.HM(5.5> 

3 . AREA, O.Dl.DH.CK.L.M, NEW, XRrIW 

FOLLOWING COMMON IS USED AS A TEMPORARY STORAGE IN THIS ROUTINE 
COMMON G(5.5) »H(5.5) 


J1=5*(M-1) 
DO 3 J=l*5 


CI=1.0 
DO 2 1=1.5 

IF ( I ,LT .5) G(If J)=CI*OF<NEW»I+l»Ul> 
H(I.J)=CK*DF(NEW» I.Jl) 


2 CI=CIM.O 

3 G ( AND*h°ARE THE COEFFICIENTS OF THE REAL AND IMAGINARY PARTS 

OF THE UPWASH 
X=DH 

DO 10 1=1. L 
JL=ML ( NEW f I ) 

IF (JL) 5.10.5 
5 Y=DH 

DO 9 J=1»JL 


Y2=Y*Y 

W(l.Jf Ilso.0 

W(2. J. I)=0»0 

YP=1.0 

DO 8 Jl=l»5 

XYP=YP 

DO 7 11=1.5 

W(1.J.I)=W(I.U.I>+G(I1.U1>*XYP 
W(2»J.I)=W(2»U»I)+H(Il#Ul) *XYP 

7 XYP=X*XYP 

8 YP=Y2*YP 

9 Y=Y+D 
10 X=X+D 

RETURN 

ENO 
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SUBROUTINE BOXPO 

C SOLUTION OF SIMULTANEOUS EQUATIONS FOR THE POTENTIAL 

COMMON A (2 *100 ,50) ,S(2+50»50) 'OA(loi5) »PS(2,5»50) »DF(2,5#50) 

1 ,ML(2' 50) , AXY (9»9)»AY(9)» XED6I (8) » YEDGl ( 8) 

2 » IE0S(2) *NS(2) »EM(50»50> , XEOG(32) »YEDG<32> ,HM(5'5> 

3 »AREA,D,DI'OH.CK»L,M»nEW,IR»IW 

C FOLLOWING COMMON IS USED AS A TEMPORARY STORAGE IN THIS ROUTINE 
COMMON E (2 '50 '51) 

11=0 

DO 9 1=1 ' L 

C ADJUST UPSTREAM INFLUENCE 
KO = 1 

NFLNS = DA (43) 

IF(NFLNS.EQ.O) GO TO 3 
KO = MAX0(1#I-NFLnS+1) 

3 continue 
JL=ML (NEW , I ) 

IF (JL.EQ.O) GO TO 9 
IF (Il.EQ.O) 60 TO 6 

C SUBTRACTION OF CONTRIBUTIONS OF PRECEDING ROWS TO UPWASH 
00 5 J=1 # JL 
DO 5 KsKO'Il 
KL=ML(NEW'K) 

K1=I+1-K 

IF (KL.EQ.O) GO To 5 
DO 4 N=1»KL 
N1=N+J 

N2=IABS(N»J)+1 

A1=A(1,N1,K1)+A(1,N2,K1) 

A2=A(2»N1,K1)+A(2,N2,K1) 

S(l'J'I)=S(l*J'I)-Al*S(l#N'K)+A2*S<2»N»K) 

4 S(2»J»I)=S(2»J'I)-A2*S(1'N»K)-A1*S(2.N»K) 

5 CONTINUE 

C SETTING UP MATRIX FOR SIMULTANEOUS EQUATIONS 

6 DO 8 J=1»UL 
DO 8 K=lrU 
N1=J+K 

N2=IABS(J-K)+1 

E(l»U»K)=A(l»Nl»l)+A(l,N2tl) 

E(2» J'K)=A(2.N1»1)+A(2.N2»1) 

E(1'K»J)=E(1'J'K) 

8 E(2»K*J)=E(2»J»K) 

C SOLUTION OF EQUATIONS 

K = MSIMEC(50'JL»l'E»S(l'l'I) ) 

IF(K.NE.I) GO TO 12 

9 11=11+1 
RETURN 

12 WRITE ( IW'41 ) 

41 FORMAT (1H010X»59HSOLUT ION OF SIMULTANEOUS EQUATIONS FOR THE POTENT 
1IAL FAILED) 

STOP 

END 
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SUBROUTINE BOXP 

COMMON A(2»100.50) ,S(2*50»50) *DA(1015) »PS(2.5»50) ,DF (2,5,50) 

1 .ML (2# 50) » AXY<9,9) .AY (9) .XEDGU8) * YEDgI < S> 

2 .IEOS(2),NS(2)*EM(50.50).XEOG(32)»YEDG(32).HM(5»5> 

3 .AREA, D.Dl.DH.CKwL.M, NEW, IR.IW 

C FOLLOWING common is used as a temporary storage In this routine 

COMMON EDG(50) .PR (2. 50) .PSI (2.50) .6(20) »X0(50) ,IXB(4) 

i .0(20.20) .8(20.2) 

DIMENSION XTEMp ( 8 ) » YTEMP ( 8 ) 

NEV=NEW 

IF (DAINEV+1002) ) 100.200*100 
100 WRITE (IW.10) 

IF(NEW.EQ.l) WRITE 
0 (Iw.12) M 

IF(NEW.EQ.2) WRITE 
0 (IW.14) M 

10 FORMAT (1H1.10X.47HUPWASH (REAL. IMAGINARY. ABSOLUTE, PHASE ANGLE)) 

12 FORMAT (1H+.58X.28H PHYSICAL PLANE MODE NO., 13/) 

14 FORMAT (1H4-.58X.31H TRANSFORMED PLANE MODE NO.»l3/) 

DO 170 1=1, L 
JL=ML(NEV.I) 

IF (JL) 110,170,110 
110 WRITE 

0 ( Iw , 20 ) I 

IF (1-1) 120,120.140 
120 DO 130 J=1 , UL 

S1=S(1.J,I) *1.57079633 
S2=S(2,J»I>*1. 57079633 
S3=S0RT < Si *Sl+S2*S2 ) 

S4=57 . 29578* AT AN2 ( S2 . S 1 ) 

130 WRITE 

0 ( Iw, 25) J,S1,S2,S3.S4 

GO TO 170 
140 JLP=JL/2 

IF(JL-2*JLP.NE.0) JLP=JLP+1 
DO 160 J=1 » JLP 

S1=SQRT(S(1»J.I)*S(1.J.I)+S(2,J,I)*S(2,J.I) ) 

S2=57 . 29578* ATAN2 ( S ( 2 , J , I ) , S ( 1 » U » I M 
J1=J+ULP 

IF(Jl.LE.JL) GO TO 150 
J1=0 
S3=0.0 
S4=0,0 
GO TO 160 
150 CONTINUE 

S3=SQRT(S<l,Ul,I)*S<l,Jl,n+S(2,Jl,I>*S(2,Ji,I)> 
S4=57.29578*ATAN2(S(2.J1»I) »S(l»Jl»I) ) 

160 WRITE 

0 (IW.25) J.S(l.J.I) ,S(2,J,D »Sl,S2»Jl,S(l»Jl»I) .S(2»J1»I) ,S3,S4 

170 CONTINUE 

C THESE ARE THE UPWASHES (ABOVE) 

200 CONTINUE 
CALL BOXPO 

C BOXPO COMPUTES THE POTENTIAL VALUES IN EACh BOX. 

C THEY ARE STORED IN THE ARRAY S. 

IF (DA(NEV+1004) ) 210,270*210 
210 WRITE (IW.30) 
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IF(NEW.EQ.l) WRITE 
0 <IW»32) M 

IF<N£W.EQ.2) WRITE 

0 (IW.34) M .. 

30 FORMAT ( 1H1 » 1 OX *6 IMPOTENT I AL CALCULATED (REAL* IMAGINARY » ABSOLUTE, 
S PHASE ANGLE)) 

32 FORMAT ( 1H+# 71X.28H PHYSICAL PLANE MODE NO. ,13/) 

34 F0RMAT(1H+»71X,31H— TRANSFORMED PLANE— MODE NO.*l3/> 

DO 250 I=1»L 
JL=ML(NEV»I> 

IF (JL) 220 » 250 # 220 
220 WRITE 

0 ( I W»20) I 

JLP=JL/2 

IF(JL-2*ULP.NE.O) JLP=ULP*1 
DO 240 U=1,JLP 

Sl=SQRT(S(l»J,I)*S(l,J#I)+S(2*J*I)*S(2»J»I) ) 
S2=57.29578*ATAN2(S<2,JtI) ,S(1#J*I> > 

J1=J+ULP 

IF(Jl.LE.JL) GO TO 230 
J1=0 
S3=0.0 
S4=0.0 
GO TO 240 
230 CONTINUE 

S3=SQRT(S(l»JI,I)*S(l,Jl*l)+S(2*>Jl»I)*S(2,Jl,I) > 

S4=57. 29578*ATAN2 ( S ( 2 , J1 # I ) , S ( 1 » Jl » I ) ) 

240 WRITE* IW, 25) J,S(l,J,I),S(2,J»I)»Sl»S2 
0 , Jl#S(l»Jl,I)*S(2*Jl# I) >S3*S4 

250 CONTINUE 

C THESE ARE THE POTENTIALS (ABOVE) 

270 CONTINUE 

C FIT OF A SERIES TO THE POTENTIAL VALUES 
300 JL=ML(NEV.L> 

NSS=NS(NEV)+l 

DY=D/YEDG(NSS) 

Y=0.5*DY 
DO 330 J=1»JL 
IF (IEoe(NEV)) 310 , 320 1 310 
310 EDG(J)= SQRT(1.0-Y*Y) 

Y=Y+DY 
GO TO 330 
320 EOG(J)=1,0 
330 CONTINUE 

N=0.5+DI*YEDG(l) 

IF (N) 360 f 360 » 340 
340 DO 350 1=1 *N 
350 X0(I)=0.0 
360 Xi=0.0 
N1=N 

Y1=YEDG(1) 

NSS=NS(NEV) 

DO 400 K=1»NSS 
X2=XEDG<K+1) 

Y2=YE0G(K+1) 

N=DI*Y2+0.5 
IF (Nl-N) 370 » 390, 390 
370 N1=N1+1 

DO 380 I=NlfN 
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Y=0*( FLOAT ( I )»0 «5) 

XO < I > =X1+ < X2-XI ) * < Y-Yl ) / < Y2-Y1 ) 

380 CONTINUE 
390 X1=X2 
Y1=Y2 
Nl-N 

400 CONTINUE 
AS=0.0 

00 430 I=1*L 
JL=ML(NEV»I) 

IF (JU 430*430*410 
410 DO 420 J=1*JL 

AS=AMAX1(AS* ABS(S(1,U*I) ) * ABSCS(2*J'I>>> 
420 CONTINUE 
430 CONTINUE 
IX=5 
IY=4 
IXY-9 

440 DO 450 1=1* IY 
450 IXB(I)=IX 
460 NC=0 

DO 470 1=1 *IY 
470 NC=NC+IXB( I ) 

DO 490 1=1 *NC 
DO 480 J=1*NC 
480 C(I*J)=0.0 
B(I*1)=0.0 
490 0(I*2)=O.O 
X1=DH 

DO 560 1=1 *L 
JL=ML(NEV,H 
Y=DH 

IF (JL) 580*580*510 
510 DO 570 J=1*JL 
XR=X1-X0(J> 

IFtXR.GE.O . ) GO TO 5l5 

IERR0R=515 

WRITE 

0 <IW*90) IERROR »NEw , NEV*I »J»X1*X0(J)»XR 

XR=0 • 5*ABS ( XR ) 

515 CONTINUE 

IF(YEDG(1) ) 520*520*530 
520 XR=XR*(XI+X0<J> J 
530 XR= SORT (XR) 

YP=1»0 

K=1 

DO 550 Nlrl.IY 

XP=XR*YP 

JX=IXB(N1) 

DO 540 N=1»JX 
G(K)=XP*EDG( J) 

XP=X1*XP 
540 K=*U1 
550 YP=Y*Y*YP 
Y=Y+D 

DO 570 N1=1*NC 
DO 560 N=1 *NC 

560 C(N1»N)=C{N1*N)+G<N1)#G(N) 


S 
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DO 570 N=l»2 

570 B(N1,N)=8(N1#N)+G(N1)*S(N* J,I) 

580 X1=X1+D 

K=MSIMER(20»NC»2*C»8) 

IF (K- 1)640 *640 *600 

C IF MSIMER FAILS# THE FOLLOWING SECTION REDUCES THE NUMBER OF TERMS 
C IN THE SERIES. 

600 DO 610 I=1*IY 
1P=I Y+l-I 

IF ( 1X8 ( IP ) ♦IP-IXY ) 610 * 620 t 620 
610 CONTINUE 
IXY=IXY-1 
GO TO 600 

620 IXB ( IP ) = IXB ( IP 1 ) -1 

IF ( IXB(IP) > 630 » 630 > 4fc0 
630 IY=IP-1 
GO TO 460 
640 K=1 

AC=0.0 

DO 670 1=1 *IY 
JX=IXB ( I ) 

CJ=1.0 

DO 660 J=1*JX 
C(K*1)=B(K*1) 

C(K*2)=B(K»2) 

AC=AMAX1(AC. ABS(C(K,1))» ABS(C(K#2))) 

IF (J.EQ.l) GO TO 660 
B ( K-l * 1 ) =C J*C ( K * 1 ) 

B ( K-l * 2 ) =C J*C ( K # 2 ) 

CJ=CJ>1.0 
660 K=K+1 

B(K-1#1)=0.0 
670 B(K-1*2>=0.0 

IF (AC-50. 0*AS) 690*690*600 1 
C GO TO 600 IF COEFFICIENTS ARE TOO LARGE. 

690 IF (DA (NEV+1006) ) 700*770*700 
C PRINTOUT OF COEFFICIENTS OF POTENTIAL 

700 CONTINUE I 

IF(NEV-l) 995 #702 .703 

702 WRITE 

0 ( IW • 42) M 

IF (NEW.EQ.l) WRITE ( Iw *45) 

IF (NEW .EQ. 2) WRITE l IW #46 ) 

GO TO 704 

703 WRITE 

0 ( IW»44) M 

704 CONTINUE 

42 FORMAT (1H1»5X»25HPHYSICAL PLANE MODE NO. #13 ) 

44 FORMAT (1H1»5X»28HTRANSF0RMED PLANE— MODE NO.* I 3/) 

45 FORMAT (1H**33X*23H (NO THICKNESS EFFECT)/) 

46 FORMAT (1H+*$3X»25H (WITH THICKNESS EFFECT)/) 

WRITE (IW*50> , 

50 FORMAT (1H-10X»53HPOTENTIAL = SUM OF P0(M»N)*X**(M*1)*Y**(2N-2)*SQR 
1TF(X) 

IF (YEDG(l)) 710»710,730 
710 WRITE (IW»52) 

IF (lEDG(NEV)) 720.750»720 
720 WRITE (IW*54) 

GO TO 750 


i 


i 
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730 WRITE < IW , 56) 

IF ( ItDG(NEV) ) 740 » 750 * 740 
740 WRITE (Iw,58) 

750 WRITE ( IW,60 ) 

52 FORMAT ( 1H+63X * 10H**2»X0**2 ) ) 

54 FORMAT ( 1H+73X * 21H*SQRTF ( 1- 1 T/YMAX) **2) > 

56 FORMAT ( 1H+63X * 4H-X0 ) ) 

58 FORMAT (1H+67X»21H*SQRTFU*<Y/YMAX)**2> ) 

60 FORMAT UH015X*52HwH£RE X = X0(Y) IS THE EQUATION 0? THE LEADING Ed 
1GE./1H020X»2IHCOEFFICIENTS P0(M,N>/1H07X» 1HM7X* 1HN14X , 9HrEAL PART 
216X*10HIMAG. PART) 

K=1 

DO 760 1=1 »IY 
JL=IXB ( I ) 

DO 760 J=1,JL 
WRITE 

0 ( IW ,65) Jfl'C(K'l) ,C(K*2) ‘ 

760 K=K+1 

WRITE (IW»70) 

770 IF(DA(NEV+1008) ) 780,860,780 
C PRINTOUT OF VALUES OF POTENTIAL AND PRESSURE 
780 WRITE (IW'80) 

IF(NEV-l) 995*782,783 

782 WRITE 

0 (IW,81) M 

IF(NEW.EQ.l) WRITE ( Iw , 83) 

IF(NEW.EQ.2) WRITE ( I W , 84) 

GO TO 784 

783 WRITE 

0 ( I W r 82 ) M 

784 WRITE(IW,85) „ . „ _ ’ 

80 F0RMAT(/1H , 10X,37HSURFACE FITTED POTENTIAL AND PRESSURE) 

81 FORMAT (1H+,47X,28H PHYSICAL PLANE MODE N0,,l3j 

82 F0RMAT(1H+,47X,31H TRANSFORMED PLANE MODE N0.,l3> 

83 FORMAT (1H+»78X,23H —NO THICKNESS EFFECT) 

84 FORMAT ( 1H+,78X,25H WITH THICKNESS EFFECT) 

85 FORMAT (1H ,22X,46H(REAL* IMAGINARY* ABSOLUTE VALUE* PHASE ANGLE)/) 
X1=DH 

DO 850 1=1, L 
«JL=ML(NEV , I ) 

Y=OH 

IF (JL) 850*850,790 
790 WRITE 

0 ( IW , 20 ) I 

DO 840 J=1,JL 
XR=X1-X0(J) 

IF(XR,GE.O.) GO TO 795 

IERR0R=795 

WRITE 

0 (IW,90) IERROR,NEw,NEV#I,J*Xl*X0(J) ,XR 

XR=0,5*ABS(XR) , 

795 CONTINUE 
XQ=0.5 

IF (YEDG(l)) 800,800,810 
800 XQ=X1 

XR=XR* (X1+X0 ( J) ) 

810 XG=XQ/XR 

XR= SORT (XR) 

DO 830 N=l*2 



PSI(N»J)=0.0 

PR(N»J)=0.0 

K=1 

YP=EDG< J) 

DO 830 N1=1»IY 

XP1=XR*YP 

JX=IXB(N1) 

DO 820 M1=1»JX 
PSI(N#J)=PSI(N,J)+C(K#N)*XPl 
PR (N» J) =PR ( N r J ) +B ( K r N ) *XP1 
XP1=X1*XP1 
820 K=K+1 
830 YP=Y*Y*YP 

PR < 1 » J) =2.0* (PR ( 1 » J) 4-PSl { t # J) *XO-PSI (2i J) *Ck ) 
PR(2 > J)=2.0*(PR(2,J)+PSI(2»J)*XQ+PSI(i.J)*CK) 
SI=SQRT(PSI(1»J)*PSI(1.J)+PSI(2^)*PSI(2»J)) 
S2=57.29578*ATAN2(PSl(2»J) #PSK1»J> ) 
S3=SQRT(PR(l»J)*PR(l,j)+PR(2»U)*PRt2»J)) 
S4=57.29578*ATAN2<PR<2,JUPR(1,U) ) 

WRIT£IIW»25) J*PSK1,J) »PSI(2t J) »S1»S2»J»PR(1»J) »PR(2»J) *S3»SU 
840 Y=Y+0 
850 Xl=Xl+D 
860 CONTINUE 

IF(NEV.EQ.l) GO TO 970 

C GET POTENTIAL FROM CORRESPONDING POINT ON TRANsFOrMAED WInG 
NEV=1 
X1=DH 

DO 950 I=1»L 
JL=ML(1»I) 

Y=DH 

IF(UL) 950*950.870 
870 DO 940 J=1»JL 
YT=Y*EM(J.I) 

JX0=YT/D ♦ 0.5 
IF( JXO.LE.O) JX0=1 
XR=X1-X0(JX0) 

XQ=0.5 

IF( YEDG(l) ) 880. 880. 890 
880 XQ=X1 

XR=XR*(X1+X0(JX0) ) 

890 CONTINUE 

IF(X1.GT.X0(JX0) ) GO TO 900 

IERROR=9QO 

WRITE 

0 ( IW. 90) IERROR,NEw,NEV»I,JtXl»X0(JX0> »XR 

XR-0 * 5*ABS ( XR ) 

900 XQ=XQ/XR 
XR=SQRT(XR) 

DO 930 N=l>2 
PSI(N*J)=0. 

K=1 

YP=EDG(JX0) 

DO 920 N1=1»IY 

XP1=XR*YP 

0X=IXB(NI) 

DO 910 Ml=l'JX 
PSI (N. J)=PSI (N» J)+C (K»N)*XPl 
XP1=X1*XP1 
910 K=K+1 
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920 yp=yt*yt*yp 

930 S(N,J*I)=PSI(N,J)/(Em<J*I)+1.E-20> 

940 Y=Y+D 
950 X1=X1+D 

00 960 1=1 #8 

IF(M.GT.1.0R.DA(26).GT.O.) GO TO 958 
XTEMP(I)=XEDG(I> 

YTEMP(I)=YEDG(I) 

958 CONTINUE 

XEDG(I)=XEDGI<I) 

960 YEDG( I )=YEDGI ( I ) 

GO TO 300 
970 CONTINUE 

IF(NEW.EQ.l) GO TO 975 

RESTORE PHYSICAL LEADING EDGES (XEDGI S YEDgI) IN XEDG $ Y^DG 
IF TRANSFORMED WING COMPUTATION FOR ALL MODES IS COMPLETED 
THIS PART OF RESTORATION HAS BEEN MADE BEFORE COMING HERE 
IF(M.GE.IFIX(DA(28>).AND.lFlX(DAt44)).LT.2) GO TO 970 
RESTORE TRANSFORMED LEADING EDGES (XTEMP * YTEmP> IN XEDG S YEDG 
TO CONTINUE OTHER MOOES COMPUTATION IN TRANSFORMED PLANE 
DO 972 1=1*8 
XEDG(I)=XTEMP(I) 

972 YEDG(I)=YTEMp(I) 

GO TO 978 

975 IF(IFIX(DA(44)).EQ.1) GO TO 978 

IF(M.LT.IFlX(DA(28)).OR.IFIX(DAt26)).EO.O> GO TO 978 

RESTORE TRANSFORMED LEAOING EDGES (XTEMP $ YTEmP) IN XEDG * YEDG 

SINCE PHYSICAL WING COMPUTATION FOR ALL MODES IS COMPLETED 

DO 977 I = 1*8 

XEDG ( I ) =XTEMP 1 1 ) 

977 YEDG ( I ) =YTEMp ( I ) 

978 CONTINUE 
JMO=5* (M-l ) 

CJ=0.0 

DO 990 1=1*5 
DO 980 U=l*5 
J0=J+dM0 
PS (1*I*JO)=0.0 
PS (2*I»JO)=0.0 
K=1 

DO 980 K2=1»IY 

J2=d+K2-1 

JX=IXB(K2) 

DO 980 Kl=l*JX 

J1=K1+I-1 

X1=AY(J2) 

IF ( J1 .GT . 1 ) X1=X1-C J*AXY ( Jl-1 * J2 ) 

Y=AXY(J1*J2)*CK 

PS (1*I*J0)=PS (1»I*J0)+ , C(K*1)*X1-C(K*2)*Y 
PS (2*I*J0)=PS (2»I*d0)+C(K(l)*Y +C(K*2»*Xl 
980 K=K+l 
990 CU=CU+1.0 
RETURN 

995 IERR0R=995 
WRITE 

0 ( IW*9Q) IERROR.NEw.NEV 

STOP 

20 FORMAT ( IH0*5X* I2*6HTH ROW I 
25 FORMAT ( 1H *5X*2(2X* 13* 1P4:13.5) > 

65 FORMAT ( 1H 218* 1P2E25.5) 
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70 FORMAT (1H1) 

90 FORMAT UH0*10X*37HBAD NUMBER IN BOXP NEAR STATEMENT N0.»l5* 
1 /1H »15X,4l5»lP5E13,5/) 

END 


SUBROUTINE MREO 

COMMON A(2»100,50) ,W<2» 50*50) »DA( 1015) »PS(2,5»50) ,DF(2,5,50> 

1 ,ML(2*50),AXY(9*9),AY(9> ,XEDGl<8)*YEDGHfl> 

2 *IEDG(2) ,NS<2> *EM(50*50) ,XEDG(32> *YEDG(32) ,HM<5*5) 

3 »AREA,D»DI*0H*CK»L,M,NEW,IR,IW 

C FOLLOWING COMMON COMMUNICATES BETWEEN *MRED* AND *LSUQA* 

COMMON XQ(150) ,YQ(i50)»DEFQ(150> »WTQ(150) »C0EC5»5) 

CONST=0. 28571429 
NP=DA(97) 

IF (NP) 83*24*30 

C A POLYNOMIAL FOR THE PrESSURE/MACH IS FITTED To VALUES 
C OF PRESSURE/MACH AT GIVEN POINTS. 

30 NX=DA(97> 

NY=DA(97) 

32 IF (100-NP) 83*35,35 
35 CONTINUE 
KP=700 

DO 11 IP=1*NP 

XQ ( IP) -DA (KP+1 ) /DA <24 } 

YQ( IP ) =DA ( KP+2 ) /DA ( 24 ) 

DEFQClP)=DA(KP+3) 

WTQ(IP)=1.0 

C DA(96)=1* INPUT DATA ARE PRESSURE COEFFICIENT 

C DA (96) =2* INPUT DATA ARE LOCAL MACH NUMBER 

IF(DA(96)-1.0> 82,106*106 

C CONVERT PRESSURE COEFFICIENT INTO LOCAL MACh NUMBER 

105 DEFQ(lP)=SQRT(5.*(1.2/( (1.+0 ,7*DEFQ ( IP) > **CONST> -1 . ) ) 

106 CONTINUE 
11 KP=KP*3 

C USE LEAST SQUARE METHOD TO CURVE FIT DATA 
IF(IFIX(DA(1011) ).NE,0) WRITE <IW»4l) 

CALL LSQUA ( NP , NX * NY * 1 0 1 1 ) 

DO 111 J=l,5 
DO 111 1=1,5 
111 HM(I,J)=COE(I*J) 

GO TO 28 

C PRESENTLY INPUT OF PRESSURE :OEFFlClENT IN 
C A POLYNOMIAL FORM IS NOT ALLOWED 

24 CONTINUE 
YP=1.0 
K=1 

DO 27 J=1 * 5 
XYP=YP 
DO 26 1=1,5 
HM(I,U)=XYP*DA(K+70) 

K=K+1 

26 XYP=XYP*DA(24) 

27 YP=YP*DA(24)**2 

28 CONTINUE 
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RETURN 

82 IPR=96 
GO TO 85 

83 IPR=97 
85 WRITE 

0 ( IW# 45) IPR 

STOP 

41 FORMAT (1H010X»56HCOMPUT£D MACH(X#Y> = SUM oF HM <N ,M) *X** (N-l ) *Y* 
l*(2M-2)/lH010X»54H(lN DIMENSIONLESS COORDINATES - DISTANCE/CHORD L 
2ENGTH) /1H09X# 1HN7X # 1HM16X# 8H HM(N#M) ) 

45 FORMAT (1H0»10X»14HMRED— BAD DATA# I 5 ) 

END 


FUNCTION CIN(X1#S) 

SINE AND COSINE INTEGRAL SUBROUTINE 
IF CALLED BY THE STATEMENT C=CIN<X»S> 

C AND S ARE THE INTEGRALS OVER T FROM 1 TO INFINITY OF 
COSUTl/T AND SIN(XT)/T 

SG=1.0 
X=X1 

IF (X) 1 »2»2 

1 SG=-SG 
X=-X 

2 X2=X*X 
IF (X-1.0) 3»3#4 

C FOR ABS(X) LESS THAN 1 A SERIES EXPANSION IS USED 

3 V=( ( (X2/98. 0-0. 6) *.05*X2+1.0)*X2/18. 0-1. 0)*X+1. 57079633 
U=( (X2/45.0-t.0)*X2/24.0+1.0)*X2/4.0-.5772l5665-ALOG(X) 

GO TO 5 

C FOR ABS(X) GREATER THAN 1 APPROXIMATIONS OF HASTINGS ARE USED 

4 P=< ( (X2+19.394119)*X2+47.41l538)*X2+8.493336)/( < < (X2+21. 361055) 

1 *X2+70. 376496) *X2+30. 038227) *X) 

Q=( ( (X2+21.383724)*X2+49.719775)*X2+5.089504)/( ( < (X2+27. 177958) 
1 *X2+U9. 918932) *X2+76#707876)*X2) 

CO=COS (X) 

SI=SIN (X) 

U=Q*CO-P*SI 

V=P*CO+Q*SI 

5 S=V*SG 
CIN=U 
RETURN 
END 
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FUNCTION MSIMERtM.NrLr AtB) 
DIMENSION A(Mrl) ) 

ERRR=1.E-14 
DO 3D I = 1>N 
c = 0.0 

DO 10 J = lrN 

io c = amaxi {c.abs(a<i.j> n 

IF (C.EQ.0.0) SO TO 1000 
00 20 J = 1*N 
20 A( I* J) = A<IrJ)/C 
DO 30 J = 1»L 
30 B ( I * J) = B(IrJ)/C 
IF(N.EQ.l) SO TO 205 
NM = N - 1 
DO 200 J = 1 rNM 
C = 0.0 
K = 0 

DO 40 I s J»N 
D = ABSlAllr J) ) 

IF (C.6E.D) 60 TO 40 
K = I 
C = D 

40 CONTINUE 

IFtK.EQ.O.OR.C.LT.ERRR) 60 TO 1000 
IF(K.EQ.J) 60 TO 70 
DO 50 JJ = JrN 
C = A l Jr JJ) 

A(J'JJ) = A IK » JJ) 

50 A(K.JJ) = C 
DO 60 JJ = lrL 
C = B( J» JJ) 

B(J*JJ) = B(KtJJ) 

60 B(Kr JJ) = C 
70 C = A(JrJ) 

JP : J M 
DO 60 J<J = JPrN 
BO A ( J • J J ) = A ( Jr JJ) /C 
90 DO 100 JJ = 1 r L 
100 B(JrJJ) = B(J»JJ)/C 
DO 200 I = 1 »N 
IF(I.EQ.J) 60 TO 200 
C = A(IrJ) 

DO UO JJ = JPrN 

110 A(IrJJ) = A(IrJJ) - C*A(JfJj) 

DO 120 JJ = 1 r L 

120 B ( I r JJ) s B(IrJJ) - C*8(J#'JJ) 

200 CONTINUE 
205 C = A(N»N) 

IF(ABS(C) .LT.ERRR) GO TO 1000 
DO 210 J s 1 r L 
210 B(N,J) = B(N* J)/C 

IF(N.EQ.l) SO TO 230 
00 220 I = 1»NM 
C = A(IrN) 

DO 220 JJ = 1»L 

220 B(IrJJ) = BUrJJ) - C*8(N*JJ) 

230 MSIMER = 1 
RETURN 
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1000 MSIMER = 2 
RETURN 
END 


FUNCTION MSIMEC(M,N.L.A»B) 

DIMENSION A(Mf 1) 

COMPLEX A. Bf G 
ERRR=1.E-14 
DO 30 I = 1*N 
C = 0.0 
DO 10 J = 1 f N 

io c=amaxi(C»abs{real(A( i ,j)i ) ,abs<aimag(a(I.J) > ) > 
if(c.eq.o.o) go to iooo 

DO 20 J = 1#N 
20 A(l'j) s A ( I f J> /C 
DO 30 J = 1 > L 
30 B ( I > J) = B ( I * J) /C 
IF(N.EQ.l) GO TO 205 
NM = N - 1 
DO 200 J S 1»NM 
C = 0.0 
K = 0 

DO 40 I = JfN 

D=ABS(REAL(A(I.J) )) +ABS< AlMAGt At I *J) ) ) 

IF(C.GE.D) GO TO 40 
K = I 
C = D 

40 CONTINUE 

IF(K.EQ.O.OR.C.LT.ERRR) GO TO 1000 
IF(K.EQ.J) GO TO 70 
DO 50 JJ = J.N 
G = A(JfJJ) 

A(J.JJ) = A(K.JJ) 

50 A(K.JJ) = G 
DO 60 JJ S l.L 
G = B(JtJJ) 

B(J»JJ) = B(KfJJ) 

60 B(K.JJ) = G 
70 G = 1<0/A( Jf J) 

JP = J ♦ 1 
DO 60 JJ = JP f N 
80 A(J.JJ) = A ( J. JJ) *G 
90 00 100 JJ = If L 
100 B(JfJJ) = Bl Jf JJ) *G 
DO 200 I = IfN 
XF(I.EQ.J) GO TO 200 
G = A ( I f J) 

DO 110 JJ = JP f N 

110 A(IfJJ) = A(IfJJ) - G*A ( J* JJ) 

DO 120 JJ = 1 f L 

120 BU.JJ) = B(ItJJ) - G*B(J«JJ> 

200 CONTINUE 
205 G = A(NfN) 

IF (ABS(REAL(G) ) ♦ ABS(AIMAG(G) ) .LT.ERRR) Go TO loOO 
DO 210 J = lfL 
210 B(NfJ) = B(NfJ)/G 
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